- Email: [email protected]
TIME DOMAIN ANALYSIS FOR DAMAGE DETECTION IN SMART STRUCTURES
Mechanical Systems and Signal Processing (1997) 11(3), 409–423
TIME DOMAIN ANALYSIS FOR DAMAGE DETECTION IN SMART STRUCTURES J. C D. J. I Center for Intelligent Material Systems and Structures, Department of Engineering Science and Mechanics, Virginia Polytechnical Institute and State University, Blacksburg, Virginia 24061-0219, U.S.A. (Received January 1995, accepted December 1996) A non-destructive time domain approach to examine structural damage is presented. Unlike frequency domain methods which generally rely on analytical models, this method is independent of modal parameters and analytical models. Time histories of the vibration response of the structure were used to identify the presence of damage. Two simple finite element models were constructed to examine axial as well as transverse vibrations. The knowledge gained from these simulations led to a series of experiments which substantiated the potential of the proposed method. Comparing time responses, e.g., displacement or velocity, due to different material defects revealed the existence of damage in cases where measured frequency shifts were minimal. 7 1997 Academic Press Limited
1. INTRODUCTION
In damage detection theory, frequency domain methods such as modal analysis have been one of the most popular approaches used to date. However, the reliability of typical experimental results in this era shown in the literature can be questioned when the damage is relatively small [1–3]. Frequency domain-related damage detection algorithms usually obtain damage characteristics by adjusting an analytical reference model to match analytical data with experimental data. As the analytical model itself can only approximate the actual mechanical structure, it introduces computational errors [3, 4]. In general, it remains difficult to detect a small change in frequency by examining transfer functions or other curve-fitting methods. However, the phenomena of beating is one that exaggerates the closeness of frequencies of two time signals. Hence, it is proposed to use the beat phenomena to look for small changes in response frequencies due to damage of complex structures. In the time domain, changes in natural frequencies due to damage result in phase shifts between the vibration response of the healthy structure and the vibration response of the damaged structure. To access this information, response functions of both damaged and undamaged structures are subtracted from one another which will in the presence of a frequency discrepancy produce a beating phenomena. Valuable results, however, can only be achieved if the experimental hardware remains unchanged throughout the entire investigation. This suggests the use of smart structures where piezoceramic patches are either surface mounted on or embedded in the structure to both actuate and measure vibrations of the host structure. The proper model of damage is itself a research topic. Damage in a structure is thought to affect the coefficients of stiffness, damping, and mass of the material(s). This work relies on a simple change of mass and stiffness as an approach to model damage. 0888–3270/97/030409 + 15 $25.00/0/mssp960086
7 1997 Academic Press Limited
410
. . .
The simulations illustrate that the time response is able to discover material damage in cases that produce very small changes in the frequency response. Note that the finite element models referred to here are not used in the damage detection algorithm. They are only used to allow simulation of vibration responses to illustrate the proposed procedure and to guide the subsequent experimental validation. By the same token, the simulations were conducted without explicit consideration of noise. Extrapolating the simulation procedure onto real structures substantiated the analytical results. For the experiments the affects of external disturbances were dismissed by utilising band-pass filters. Specific noise analysis is not part of this research. Moving from a plain aluminum plate to a section of a helicopter rotor blade illustrates that the procedure works on complex structures as well. 2. BEAT METHOD
The damage detection method presented here involves excitation of the structure with the time history mentioned above and measurement of the time response at various intervals, much like a monitoring approach. The system is excited until steady state is reached, then released to let it vibrate in free decay. This motion is dominated by a specific mode depending on the input signal. To select the changes of natural frequencies due to damage without reference to the analytical model or the frequency response, the behaviour of the healthy structure is investigated to provide the reference data. It is assumed that the input signal as well as the location of sensors and actuators remain unchanged during the entire test (which for the later experiments was achieved by using surface-bonded piezoceramic patches). If there are any material defects, the corresponding time responses from before and after the damage should then differ by a certain phase shift. The phase shift itself reveals the changes in the dominant natural frequency which are caused through damage. A simplified but effective method to identify this shift in phase, is to subtract the corresponding time responses from each other shown in Fig. 1. As long as the frequencies do not differ by more than approximately 30%, the result of this subtraction will be a beat phenomena. In fact, the time for each beat is the time the responses need to reach a 2p phase shift. The less significant the damage, the smaller the frequency changes and the lower the beating frequency will be. As expected, for the undamaged case the time for a 2p phase shift, t360 , is infinite. The quantity t360 determines how critical the damage is to the structure’s response, and is the focus of this work. This method can be very time consuming, especially in low
Figure 1. Time domain detection algorithm.
411
Figure 2. Multifrequency burst input signal.
damage cases. To improve the time factor, the quantity t180 is defined. It is derived by subtracting the absolute values of the corresponding time responses. Du(t) = =uh (t)= − =ud (t)=
(1)
The phase shift is exactly 180° when this difference becomes zero. It is assumed that both time functions vibrate in only one mode with a similar amplitude. For small enough damage, equation (1) then results in beating with periodic minimums every t180 s. Equation (1) therefore defines the detection algorithm proposed here. To reveal a frequency shift between two time functions from subtraction, each function must be dominated by only one frequency. The signal shown in Fig. 2 is a combination of both impulse signal and sinusoidal waves. It allows excitation of any designated frequency band. Preparing this signal to cover a frequency range around the mode of interest reduces the magnitude’s dependence on changes in the natural frequency due to damage. The system will vibrate approximately with its maximum amplitude, which is also unaffected by the time the excitation is stopped. This input signal is called Schroeder signal [5].
3. NUMERICAL VERIFICATION
3.1. Different degrees of damage were simulated and detected by exciting an analytical finite element model with the Schroeder input signal of Fig. 2. The natural frequencies for axial vibration referring to the continuous fixed-free model are known to be fn =
(2n − 1) 2l
X
E with n = 1, 2, . . . r
The properties of the bar model to be examined were taken to be A = 10 mm2 l = 1000 mm
(2)
412
. . .
Figure 3. First mode time domain analysis for a 1% mass and stiffness loss at element 3 of the bar model.
E = 7.0 × 107
kg mm s2
r = 2.7 × 10−6
kg mm3
(3)
which then causes the frequencies of equation (2) to be fn = (2n − 1) × 1273 Hz with n = 1, 2, . . .
(4)
Using a 10-element model (Appendices A and B) provided a reasonable convergence to the continuous model, since the dominant natural frequencies correspond fairly well, i.e. f1 = 1274 Hz + 0.07% f2 = 3854 Hz + 0.90%
(5)
Based on equation (B6) the first and the second mode of the system are then excited with the Schroeder signal for 12.5 ms. Figures 3 and 4 show the beating functions obtained by the detection algorithm explained in Section 2. The frequency bandwidth of the applied input in both cases covers a total range of 50 Hz, 25 Hz above and 25 Hz below the respective natural frequency of
Figure 4. Second mode time domain analysis for a 1% mass and stiffness loss at element 3 of the bar model.
413
the undamaged bar model. The damage assigned reduces the mass and stiffness of element three, counted from the clamp, of the bar model by only 1%. The total mass and stiffness loss for the 10-element model would then be 0.1%. The system was designed with very little damping using equations (B2–B5); the damping ratios for the first and second natural frequency are 1.25 per thousand. Both modes examined above are affected by the simulated stiffness and mass reduction in a similar fashion. Each natural frequency shifts to a slightly different value when damage occurs. The magnitude of the beating function representing the second mode is approximately one-third the value of the first mode analysis. Since the actual reading is the time for a single beat, the magnitude in this respect is insignificant and therefore was adjusted to allow optimal measurement of t180 . 3.2. The simulation with the beam model derived in Appendix A differs slightly from the bar model analysis. The structural properties of the specimen [equation (3)] remain unchanged, yet the bending motion produces natural frequencies significantly smaller than the axial motion. The natural frequencies for the continuous model of a clamped-free beam in transverse vibration are known as
fn =
0 1X
1 kn 2p l
2
EI with n = 1, 2, . . . rA
(6)
whereas kn is the weighted natural frequency bn l is taken from [6]. Inserting values into equation (6) leads to
fn = 2.3394 × kn2 with kn =
8
bn l
nE5
(2n − 1)p n q 5. 2
(7)
The second moment of area is taken from equation (A19), with b = 1 mm and h = 10 mm. The 10-element model converged well with the continuous model; the first and second natural frequencies were computed to be f1 = 8.21 Hz + 0.22% f2 = 51.57 Hz + 0.04%
(8)
As before the dyanamic model defined in equation (B6) was used to simulate the system’s time responses. For the bending motion, the damping ratio of the first two modes was set to 12.5 per thousand. This time the frequency bandwidth of the Schroeder input covered approximately a range of 5 Hz, symmetric to the respective natural frequency of interest; the structure was excited for 0.5 s. The damage assigned, reduces the stiffness of a designated element per definition and calculates through equation (A20) the corresponding loss of mass for the same element. Figures 5 and 6 give examples of beating functions obtained from this simulation. Again t180 is easy to measure and reveals clearly the presence of damage.
414
. . .
Figure 5. First mode time domain analysis for a 1% stiffness and a 0.34% mass reduction at element 3 of the bar model.
4. EXPERIMENTAL VERIFICATION
4.1. Figure 7 shows a schematic of the cantilevered aluminum plate with piezoceramic sensors and actuators bonded to the plate’s surface for excitation and measurement. The damage is applied in the form of a concentrated mass that is fastened to the upper edge of the specimen. Three different Schroeder-phased input signals to excite the three lowest modes of the structure separately were generated through the Matlab Simulink Toolbox. The signals were then transferred via DSpace computer to a Tektronix Analyzer 2600 unit which allowed the data to be stored in a PC. Sending the recorded data as output through the analyser to excite the plate guaranteed repeatable excitation. Each excitation signal contained 512 data points and of the vibration response a total of 4096 data points were recorded after the excitation was shut off. The sampling rate of the analyser was set to approximately 2.5 times the frequency of the active mode. To tackle the undesired effect of noise, the response signals were applied to a series of ITHACO 4302 Dual 24 dB/Octave low- and high-pass filters. In conjunction, they produced the effect of band-pass filtering. For each set of excitation signals, the filters were adjusted as to optimise the response in the desired frequency range and to dismiss unwanted noise. In the course of the experiments, no additional investigation of the
Figure 6. Second mode time domain analysis for a 1% stiffness and a 0.34% mass reduction at element 6 of the bar model.
415
Figure 7. Plate, experimental set-up.
affect of induced noise was conducted. The first mode of the plate with a frequency of approximately 3 Hz was excited for 40 s with a Schroeder-phased signal covering a frequency range of 2.5–3.5 Hz. The resulting free decay was recorded for 80 s, once for the ‘undamaged’ case as a reference signal, and several times with the mass at different positions along the upper edge of the plate. Figure 8 shows how the mass at position 10 affected the first mode. It forced the frequency of this mode to shift a minimal 0.5%. This information can be retrieved from the size of the beats. The time extent of each beat is inversely proportional to the shift in phase where the factor of proportional depends on the corresponding natural frequency: Df = c ·
1 1 , c = , i = mode 2fi t180
(9)
To examine the second mode, f2 3 9 Hz, the structure was excited for 10 s with a 7.5–9.5-Hz Schroeder-phase input signal and the response was recorded for 16 s: again, running one ‘undamaged’ case first and then several cases with the mass attached at
Figure 8. First mode response analysis.
416
. . .
Figure 9. Second mode response analysis.
different positions. With the mass attached at position 5, closer to the clamp than the previous case, the measured time for a 180° phase shift was 3.5 s which led to a change in frequency of 1.6% (Fig. 9). Mode three, f3 3 18.5 Hz, provided similar results, but it also revealed how damping can interfere with retrieving the desired data. The input signal was applied for 10 s covering a frequency range of 17–20 Hz; the response was recorded for 16 s. Again, the mass was attached at location 5 (Fig. 10). Figures 9 and 10 display how modes 2 and 3 were affected by certain changes of the structure. For this particular case, the first mode’s time response measurements did not reveal a frequency shift in time. This leads to the conclusion that this particular eigenfrequency was not affected at all by the damage. This series of measurements can be very valuable to damage localisation procedures which consider rank-ordering of shifts in eigenfrequencies.
Figure 10. Third mode response analysis.
417
Embedded PZT
Figure 11. Rotor blade section.
4.2. The section of a helicopter rotor blade investigated in this section measures 0.546 × 0.533 m with a maximum thickness of 0.05 m. The centre core of the structure, an oval hollow shaft with a wall thickness of 3.2 mm, is a metal alloy giving support to the blade profile which is made of composite materials. The blade was clamped to a concrete monolith (Fig. 11). Two piezoceramic devices each 9.5 mm square were glued to the blade surface above the metal core, one for excitation and one for measurement. To introduce the non-destructive damage (thus allowing repeatability of experiments), a concentrated mass was attached to the surface at different locations. The mass added in this experiment was about 0.1% of the total mass of the structure. To find frequencies of resonance, the structure was exposed to chirp signals covering certain frequency bands and the intensity of the response was monitored. The two lowest frequencies of strongest response were found at 1.25 and 1.8 kHz. Two narrowband Schroeder-phased signals were generated for each of these two frequencies and stored to a PC following the procedure outlined for the previous experiment. The recorded signals were sent through the Tektronix Analyzer to excite the structure for 150 ms and the free decay response was monitored for 850 ms sampling 4096 data points. Figures 12 and 13 show examples of the detected beating phenomena and demonstrate how both frequencies are affected by the local increase of mass. Not referring to frequency domain functions, it is unknown which modes are actually being measured. However, this knowledge is expendable in the presence of damage and it is of no real advantage to this damage detection approach. In Fig. 12 the reference function oscillates at a frequency of
Figure 12. Blade response analysis.
418
. . .
Figure 13. Blade response analysis.
1250 Hz and the added mass causes a phase shift of 7.7%. Here, the measure of one beat becomes critical due to the high damage and the limited memory availability of the storage unit used in this experiment. The reference frequency in Fig. 13 was measured to be 1800 Hz and the effect of the added mass forced a shift in phase of 15%. Whereas Figs 12 and 13 are an advantageous way to visualise the presence of damage, Fig. 14 gives more insight to where a shift in phase actually occurs and allows for more exact measurement of the described beats of Fig. 13. The altered frequency takes exactly 18 ms to outrun the reference signal by half a period. Also noticeable is the decrease in amplitude due to the damping of the composite structure. 5. CONCLUSIONS
A time domain procedure has been presented which relies on the measured time response due to a simple harmonic input of a healthy, yet unmodeled structure. The healthy signal is then compared to the measured time response for identical input of the structure after some local damage is thought to have occurred. The comparison is attained by subtracting the two signals from one another. The resulting beat phenomena provides an indication of the existence and extent of damage reflected in local mass and/or stiffness changes.
Figure 14. Third mode response analysis. – – – –, Damaged; –––, healthy.
419
Numerical simulations and tests on a plate-like structure and an actual helicopter blade section have illustrated the algorithms capability of identifying changes in frequencies as small as 0.5%. A small reduction (hole) or increase (bullet) in mass can be detected experimentally with a minimum of difficulty and absolutely no a priori information about the structure. Damage detection algorithms without a priori knowledge of any model and without first estimating a model are potentially useful for informing maintenance personnel of the existence of damage in a variety of civilian and military applications. These applications range from civil structures and infrastructures to manufacturing product imperfections to helicopter components. The simplicity of the algorithm allows on board diagnostics when coupled with self-sensing piezoceramic embedded composites. ACKNOWLEDGEMENTS
The authors gratefully acknowledge the support of ARO grant number DAAL 03-92-G-0180 under the supervision of Dr Gary Anderson. REFERENCES 1. P. F. R and N. A 1990 Journal of Sound and Vibration 138, 381–388. Identification of crack location and magnitude in a cantilevered beam from the vibration modes. 2. C. D, P. Q. Z, W. Q. F and T. C. H 1994 Proceedings of the 12th IMAC, Honolulu, HW, 98–105. The sensitivity study of the modal parameters of a cracked beam. 3. H. T. B, D. J. I, D. J. L and Y. W 1996 Journal of Sound and Vibration 191, 859–880. An experimentally validated damage detection theory in smart structures. 4. J. C, G. R. T and K. W 1994 Proceedings of the 12th IMAC, Honolulu, HW, 778–786. A simplified approach to the numerical and experimental modelling of the dynamics of a cracked beam. 5. Y. Y, D. S. B and R. E. S 1993 ACC Proceedings, San Francisco, CA, 3021–3023. Frequency domain identification for robust large space structures control design. 6. D. J. I 1994 Engineering Vibration. Englewood Cliffs, NJ: Simon and Schuster. 7. D. K. L and G. K 1994 AIAA Proceedings, Hilton Head, SC, 192–198. Location and estimation of damage in a beam using different identification algorithms.
APPENDIX A
A first step in assessing the technique’s feasibility is to set up a simulation for the longitudinal vibration of a cantilevered bar. The bar is divided into ne elements of length si allowing each node one single degree of freedom [6]. The stiffness matrix for the ith element is Ki =
$
Ei Ai 1 si −1
%
−1 EA = i i K i 1 si
(A1)
where Ei is Young’s modulus and Ai the cross-sectional area. The mass matrix is Mi =
$ %
rAi si 2 6 1
1 rA s = i i M i 2 6
(A2)
where r is the mass density. Since each element is attached at the ith and (i + 1)th nodes,
. . .
420
the global matrices can be obtained by direct assembly. Considering the algorithms provided in [7] the expanded stiffness matrix in global co-ordinates is Khi =
Ei Ai Bi K i BiT si
(A3)
introducing the matrix Bi Bi = [ei ei + 1 ]
(A4)
where ej is the jth unit vector, i.e. all elements are zero except for a value 1 in the jth position. Summing the individual stiffness matrices results in ne
ne Ei Ai Bi K i BiT = s Khi . si i=1 i=1
Kh = s
(A5)
The global mass matrix is assembled in exactly the same manner with rAi si Bi M i BiT 6
(A6)
ne rAi si Bi M i BiT = s Mhi . 6 i=1 i=1
(A7)
Mhi = which yields ne
Mh = s
The subscript h denotes quantities from the healthy model. The damage is represented by changes of the cross-sectional area which equal losses of stiffness and mass. While this is not the most advanced model of damage, it does provide a reasonable set of data in a controlled environment for use in testing the approach. Note that such a model is just a subset of all damage models. The cross-section for the ith element Ai will be varied by the scalar factor ai (1 − ai )Ai for 0 E ai E 1
(A8)
which is called the damage coefficient of the ith element of the beam. Suppose that any element may have sustained damage as sketched in Fig. A1. The damaged global stiffness and mass matrices then will change to ne
ne Ei (1 − ai )Ai Bi K i BiT = s (1 − ai )Kdi si i=1 i=1
Kd = s
Figure A1. Bar with reduction in cross-sectional area.
(A9)
ne ne r(1 − ai )Ai si Md = s Bi M i BiT = s (1 − ai )Mdi 6 i=1 i=1
421 (A10)
The subscript d denotes quantities from the damaged model. The beam model was set up to develop further the results obtained by axial excitations. Referring to Euler–Bernoulli beam elements, stiffness and mass matrices are given by equations (A11) and (A12). All constants are defined as before and Ii denotes the second moment of inertia. 6si −12 6si L K 12 2 G 6si 4si −6si 2si2 G Ei Ii Ki = 3 G G = Esi3Ii K i si −12 −6si 12 −6si i G G k 6si 2si2 −6si 4si2 l
(A11)
22si 54 −13si L K 156 G 22si 4si2 13si −3si2 G rA s i si Mi = i iG G = rA M . 420 54 13si 156 −22si 420 i G G k −13si −3si2 −22si 4si2 l
(A12)
Just as in Section 2.1, the global stiffness and global mass matrices are given by equations (A5) and (A7), respectively. The transformation matrix appears a little different from (A4) due to two degrees of freedom at each node. In particular Bi = [e2i − 1 e2i
e2i + 1
e2i + 2 ].
(A13)
The damage in a beam element is defined as a change in the flexural rigidity Ei Ii that will not affect the mass of the beam. It is assumed also that only sustained damage is modelled. Following the example of the bar, a damage coefficient decreases the values of the corresponding matrix elements. In this case, only the stiffness matrix is altered as indicated (1 − di )
Ei Ii K for 0 E di E 1. si3 i
(A14)
The global stiffness matrix then becomes ne
ne (1 − di )Ei Ii Bi K i BiT = s (1 − di )Kdi 3 si i=1 i=1
Kd = s
(A15)
and the global mass matrix remains ne ne rAi si Bi M i BiT = s Mdi = s Mhi . 420 i=1 i=1 i=1 ne
Md = Mh = s
(A16)
Suppose the damage is modeled as a crack such as a cut with depth Dhi shown in Fig. A2. In that case, it is obvious that of all constants involved in equation (A11) the
. . .
422
Figure A2. Bar with reduction in cross-sectional area.
damage coefficient di to represent the crack it is possible to derive a correlation between the depth of the crack and the percentage of stiffness reduction. Thus Ixx (Dhi ) = (1 − di )Ixx .
(A17)
The second moment of the area of the cracked element about the x-axis is denoted as Ixx (Dhi ) =
g
h/2 − Dhi
y 2bdy
−h/2
= 121 bh 3 − 14 bh 2Dhi + 12 bhDhi2 − 13 bDhi3
(A18)
and of the area of the healthy elements Ixx = 121 bh 3,
(A19)
respectively. Using equations (A18) and (A19) in (A17) results in
0 1 0 1 0 1 Dhi h
3
− 32
Dhi h
2
+ 34
Dhi − 14 di = 0. h
(A20)
Since the ratio Dhi /h is proportional to the loss of mass of the cracked element, equation (A20) evaluates the correlation of stiffness reduction and mass reduction for this one particular element. In Fig. A3 the relation between the two damage ratios di and Dhi /h is plotted. The assumption that damage, as modeled before, will not alter the mass but stiffness exclusively, needs to be applied carefully. It is shown that there is an
Figure A3. Damage ratio correlation.
423
immediate interrelation between the two quantities when damage is modeled as suggested in Fig. A2. Due to this result the global mass matrix defined in equation (A16) will now become ne
Md = s i=1
ne
=s i=1
0
1−
Dhi rAi si B M B T h 420 i i i
1
0
1−
Dhi Mdi with Dhi = f(di ) h
1
(A21)
APPENDIX B
The dynamic model for both undamaged mechanical structures, beam and bar is given by Mh u¨ (t) + Cu˙ (t) + Kh u(t) = F(t).
(B1)
The damping given in equation (B1) is assumed to be proportional, depending on both global stiffness and global mass matrices of the form C = aMh + bKh .
(B2)
Here, the constants a and b may be defined freely, yet it remains difficult to understand how the proportional damping will actually affect the model’s response. The damping ratio zj =
a bv + j 2vj 2
(B3)
takes remedial action. It determines the rate of decay of the jth natural frequency. Using expression (B3) derived in [6] in equation (B2) leads to a = 2v1 v2 b=2
z1 v2 − z2 v1 v22 − v12
z2 v2 − z1 v1 v22 − v12
(B4) (B5)
where the first and second natural frequencies are chosen, in regard to their significant influence on the system’s response. Substituting equation (B4) and (B5) into (B2) provides a more transparent damping behaviour. In the case of damage the dynamic behaviour will be modeled by Md u¨ (t) + Cu˙ (t) + Kd u(t) = F(t).
(B6)
Damping will be calculated as for the undamaged system but using Kd and Md , the global matrices of the damaged model. Note again that these models are needed to provide data and are not part of the ‘detection’ algorithm proposed here.
TIME DOMAIN ANALYSIS FOR DAMAGE DETECTION IN SMART STRUCTURES J. C D. J. I Center for Intelligent Material Systems and Structures, Department of Engineering Science and Mechanics, Virginia Polytechnical Institute and State University, Blacksburg, Virginia 24061-0219, U.S.A. (Received January 1995, accepted December 1996) A non-destructive time domain approach to examine structural damage is presented. Unlike frequency domain methods which generally rely on analytical models, this method is independent of modal parameters and analytical models. Time histories of the vibration response of the structure were used to identify the presence of damage. Two simple finite element models were constructed to examine axial as well as transverse vibrations. The knowledge gained from these simulations led to a series of experiments which substantiated the potential of the proposed method. Comparing time responses, e.g., displacement or velocity, due to different material defects revealed the existence of damage in cases where measured frequency shifts were minimal. 7 1997 Academic Press Limited
1. INTRODUCTION
In damage detection theory, frequency domain methods such as modal analysis have been one of the most popular approaches used to date. However, the reliability of typical experimental results in this era shown in the literature can be questioned when the damage is relatively small [1–3]. Frequency domain-related damage detection algorithms usually obtain damage characteristics by adjusting an analytical reference model to match analytical data with experimental data. As the analytical model itself can only approximate the actual mechanical structure, it introduces computational errors [3, 4]. In general, it remains difficult to detect a small change in frequency by examining transfer functions or other curve-fitting methods. However, the phenomena of beating is one that exaggerates the closeness of frequencies of two time signals. Hence, it is proposed to use the beat phenomena to look for small changes in response frequencies due to damage of complex structures. In the time domain, changes in natural frequencies due to damage result in phase shifts between the vibration response of the healthy structure and the vibration response of the damaged structure. To access this information, response functions of both damaged and undamaged structures are subtracted from one another which will in the presence of a frequency discrepancy produce a beating phenomena. Valuable results, however, can only be achieved if the experimental hardware remains unchanged throughout the entire investigation. This suggests the use of smart structures where piezoceramic patches are either surface mounted on or embedded in the structure to both actuate and measure vibrations of the host structure. The proper model of damage is itself a research topic. Damage in a structure is thought to affect the coefficients of stiffness, damping, and mass of the material(s). This work relies on a simple change of mass and stiffness as an approach to model damage. 0888–3270/97/030409 + 15 $25.00/0/mssp960086
7 1997 Academic Press Limited
410
. . .
The simulations illustrate that the time response is able to discover material damage in cases that produce very small changes in the frequency response. Note that the finite element models referred to here are not used in the damage detection algorithm. They are only used to allow simulation of vibration responses to illustrate the proposed procedure and to guide the subsequent experimental validation. By the same token, the simulations were conducted without explicit consideration of noise. Extrapolating the simulation procedure onto real structures substantiated the analytical results. For the experiments the affects of external disturbances were dismissed by utilising band-pass filters. Specific noise analysis is not part of this research. Moving from a plain aluminum plate to a section of a helicopter rotor blade illustrates that the procedure works on complex structures as well. 2. BEAT METHOD
The damage detection method presented here involves excitation of the structure with the time history mentioned above and measurement of the time response at various intervals, much like a monitoring approach. The system is excited until steady state is reached, then released to let it vibrate in free decay. This motion is dominated by a specific mode depending on the input signal. To select the changes of natural frequencies due to damage without reference to the analytical model or the frequency response, the behaviour of the healthy structure is investigated to provide the reference data. It is assumed that the input signal as well as the location of sensors and actuators remain unchanged during the entire test (which for the later experiments was achieved by using surface-bonded piezoceramic patches). If there are any material defects, the corresponding time responses from before and after the damage should then differ by a certain phase shift. The phase shift itself reveals the changes in the dominant natural frequency which are caused through damage. A simplified but effective method to identify this shift in phase, is to subtract the corresponding time responses from each other shown in Fig. 1. As long as the frequencies do not differ by more than approximately 30%, the result of this subtraction will be a beat phenomena. In fact, the time for each beat is the time the responses need to reach a 2p phase shift. The less significant the damage, the smaller the frequency changes and the lower the beating frequency will be. As expected, for the undamaged case the time for a 2p phase shift, t360 , is infinite. The quantity t360 determines how critical the damage is to the structure’s response, and is the focus of this work. This method can be very time consuming, especially in low
Figure 1. Time domain detection algorithm.
411
Figure 2. Multifrequency burst input signal.
damage cases. To improve the time factor, the quantity t180 is defined. It is derived by subtracting the absolute values of the corresponding time responses. Du(t) = =uh (t)= − =ud (t)=
(1)
The phase shift is exactly 180° when this difference becomes zero. It is assumed that both time functions vibrate in only one mode with a similar amplitude. For small enough damage, equation (1) then results in beating with periodic minimums every t180 s. Equation (1) therefore defines the detection algorithm proposed here. To reveal a frequency shift between two time functions from subtraction, each function must be dominated by only one frequency. The signal shown in Fig. 2 is a combination of both impulse signal and sinusoidal waves. It allows excitation of any designated frequency band. Preparing this signal to cover a frequency range around the mode of interest reduces the magnitude’s dependence on changes in the natural frequency due to damage. The system will vibrate approximately with its maximum amplitude, which is also unaffected by the time the excitation is stopped. This input signal is called Schroeder signal [5].
3. NUMERICAL VERIFICATION
3.1. Different degrees of damage were simulated and detected by exciting an analytical finite element model with the Schroeder input signal of Fig. 2. The natural frequencies for axial vibration referring to the continuous fixed-free model are known to be fn =
(2n − 1) 2l
X
E with n = 1, 2, . . . r
The properties of the bar model to be examined were taken to be A = 10 mm2 l = 1000 mm
(2)
412
. . .
Figure 3. First mode time domain analysis for a 1% mass and stiffness loss at element 3 of the bar model.
E = 7.0 × 107
kg mm s2
r = 2.7 × 10−6
kg mm3
(3)
which then causes the frequencies of equation (2) to be fn = (2n − 1) × 1273 Hz with n = 1, 2, . . .
(4)
Using a 10-element model (Appendices A and B) provided a reasonable convergence to the continuous model, since the dominant natural frequencies correspond fairly well, i.e. f1 = 1274 Hz + 0.07% f2 = 3854 Hz + 0.90%
(5)
Based on equation (B6) the first and the second mode of the system are then excited with the Schroeder signal for 12.5 ms. Figures 3 and 4 show the beating functions obtained by the detection algorithm explained in Section 2. The frequency bandwidth of the applied input in both cases covers a total range of 50 Hz, 25 Hz above and 25 Hz below the respective natural frequency of
Figure 4. Second mode time domain analysis for a 1% mass and stiffness loss at element 3 of the bar model.
413
the undamaged bar model. The damage assigned reduces the mass and stiffness of element three, counted from the clamp, of the bar model by only 1%. The total mass and stiffness loss for the 10-element model would then be 0.1%. The system was designed with very little damping using equations (B2–B5); the damping ratios for the first and second natural frequency are 1.25 per thousand. Both modes examined above are affected by the simulated stiffness and mass reduction in a similar fashion. Each natural frequency shifts to a slightly different value when damage occurs. The magnitude of the beating function representing the second mode is approximately one-third the value of the first mode analysis. Since the actual reading is the time for a single beat, the magnitude in this respect is insignificant and therefore was adjusted to allow optimal measurement of t180 . 3.2. The simulation with the beam model derived in Appendix A differs slightly from the bar model analysis. The structural properties of the specimen [equation (3)] remain unchanged, yet the bending motion produces natural frequencies significantly smaller than the axial motion. The natural frequencies for the continuous model of a clamped-free beam in transverse vibration are known as
fn =
0 1X
1 kn 2p l
2
EI with n = 1, 2, . . . rA
(6)
whereas kn is the weighted natural frequency bn l is taken from [6]. Inserting values into equation (6) leads to
fn = 2.3394 × kn2 with kn =
8
bn l
nE5
(2n − 1)p n q 5. 2
(7)
The second moment of area is taken from equation (A19), with b = 1 mm and h = 10 mm. The 10-element model converged well with the continuous model; the first and second natural frequencies were computed to be f1 = 8.21 Hz + 0.22% f2 = 51.57 Hz + 0.04%
(8)
As before the dyanamic model defined in equation (B6) was used to simulate the system’s time responses. For the bending motion, the damping ratio of the first two modes was set to 12.5 per thousand. This time the frequency bandwidth of the Schroeder input covered approximately a range of 5 Hz, symmetric to the respective natural frequency of interest; the structure was excited for 0.5 s. The damage assigned, reduces the stiffness of a designated element per definition and calculates through equation (A20) the corresponding loss of mass for the same element. Figures 5 and 6 give examples of beating functions obtained from this simulation. Again t180 is easy to measure and reveals clearly the presence of damage.
414
. . .
Figure 5. First mode time domain analysis for a 1% stiffness and a 0.34% mass reduction at element 3 of the bar model.
4. EXPERIMENTAL VERIFICATION
4.1. Figure 7 shows a schematic of the cantilevered aluminum plate with piezoceramic sensors and actuators bonded to the plate’s surface for excitation and measurement. The damage is applied in the form of a concentrated mass that is fastened to the upper edge of the specimen. Three different Schroeder-phased input signals to excite the three lowest modes of the structure separately were generated through the Matlab Simulink Toolbox. The signals were then transferred via DSpace computer to a Tektronix Analyzer 2600 unit which allowed the data to be stored in a PC. Sending the recorded data as output through the analyser to excite the plate guaranteed repeatable excitation. Each excitation signal contained 512 data points and of the vibration response a total of 4096 data points were recorded after the excitation was shut off. The sampling rate of the analyser was set to approximately 2.5 times the frequency of the active mode. To tackle the undesired effect of noise, the response signals were applied to a series of ITHACO 4302 Dual 24 dB/Octave low- and high-pass filters. In conjunction, they produced the effect of band-pass filtering. For each set of excitation signals, the filters were adjusted as to optimise the response in the desired frequency range and to dismiss unwanted noise. In the course of the experiments, no additional investigation of the
Figure 6. Second mode time domain analysis for a 1% stiffness and a 0.34% mass reduction at element 6 of the bar model.
415
Figure 7. Plate, experimental set-up.
affect of induced noise was conducted. The first mode of the plate with a frequency of approximately 3 Hz was excited for 40 s with a Schroeder-phased signal covering a frequency range of 2.5–3.5 Hz. The resulting free decay was recorded for 80 s, once for the ‘undamaged’ case as a reference signal, and several times with the mass at different positions along the upper edge of the plate. Figure 8 shows how the mass at position 10 affected the first mode. It forced the frequency of this mode to shift a minimal 0.5%. This information can be retrieved from the size of the beats. The time extent of each beat is inversely proportional to the shift in phase where the factor of proportional depends on the corresponding natural frequency: Df = c ·
1 1 , c = , i = mode 2fi t180
(9)
To examine the second mode, f2 3 9 Hz, the structure was excited for 10 s with a 7.5–9.5-Hz Schroeder-phase input signal and the response was recorded for 16 s: again, running one ‘undamaged’ case first and then several cases with the mass attached at
Figure 8. First mode response analysis.
416
. . .
Figure 9. Second mode response analysis.
different positions. With the mass attached at position 5, closer to the clamp than the previous case, the measured time for a 180° phase shift was 3.5 s which led to a change in frequency of 1.6% (Fig. 9). Mode three, f3 3 18.5 Hz, provided similar results, but it also revealed how damping can interfere with retrieving the desired data. The input signal was applied for 10 s covering a frequency range of 17–20 Hz; the response was recorded for 16 s. Again, the mass was attached at location 5 (Fig. 10). Figures 9 and 10 display how modes 2 and 3 were affected by certain changes of the structure. For this particular case, the first mode’s time response measurements did not reveal a frequency shift in time. This leads to the conclusion that this particular eigenfrequency was not affected at all by the damage. This series of measurements can be very valuable to damage localisation procedures which consider rank-ordering of shifts in eigenfrequencies.
Figure 10. Third mode response analysis.
417
Embedded PZT
Figure 11. Rotor blade section.
4.2. The section of a helicopter rotor blade investigated in this section measures 0.546 × 0.533 m with a maximum thickness of 0.05 m. The centre core of the structure, an oval hollow shaft with a wall thickness of 3.2 mm, is a metal alloy giving support to the blade profile which is made of composite materials. The blade was clamped to a concrete monolith (Fig. 11). Two piezoceramic devices each 9.5 mm square were glued to the blade surface above the metal core, one for excitation and one for measurement. To introduce the non-destructive damage (thus allowing repeatability of experiments), a concentrated mass was attached to the surface at different locations. The mass added in this experiment was about 0.1% of the total mass of the structure. To find frequencies of resonance, the structure was exposed to chirp signals covering certain frequency bands and the intensity of the response was monitored. The two lowest frequencies of strongest response were found at 1.25 and 1.8 kHz. Two narrowband Schroeder-phased signals were generated for each of these two frequencies and stored to a PC following the procedure outlined for the previous experiment. The recorded signals were sent through the Tektronix Analyzer to excite the structure for 150 ms and the free decay response was monitored for 850 ms sampling 4096 data points. Figures 12 and 13 show examples of the detected beating phenomena and demonstrate how both frequencies are affected by the local increase of mass. Not referring to frequency domain functions, it is unknown which modes are actually being measured. However, this knowledge is expendable in the presence of damage and it is of no real advantage to this damage detection approach. In Fig. 12 the reference function oscillates at a frequency of
Figure 12. Blade response analysis.
418
. . .
Figure 13. Blade response analysis.
1250 Hz and the added mass causes a phase shift of 7.7%. Here, the measure of one beat becomes critical due to the high damage and the limited memory availability of the storage unit used in this experiment. The reference frequency in Fig. 13 was measured to be 1800 Hz and the effect of the added mass forced a shift in phase of 15%. Whereas Figs 12 and 13 are an advantageous way to visualise the presence of damage, Fig. 14 gives more insight to where a shift in phase actually occurs and allows for more exact measurement of the described beats of Fig. 13. The altered frequency takes exactly 18 ms to outrun the reference signal by half a period. Also noticeable is the decrease in amplitude due to the damping of the composite structure. 5. CONCLUSIONS
A time domain procedure has been presented which relies on the measured time response due to a simple harmonic input of a healthy, yet unmodeled structure. The healthy signal is then compared to the measured time response for identical input of the structure after some local damage is thought to have occurred. The comparison is attained by subtracting the two signals from one another. The resulting beat phenomena provides an indication of the existence and extent of damage reflected in local mass and/or stiffness changes.
Figure 14. Third mode response analysis. – – – –, Damaged; –––, healthy.
419
Numerical simulations and tests on a plate-like structure and an actual helicopter blade section have illustrated the algorithms capability of identifying changes in frequencies as small as 0.5%. A small reduction (hole) or increase (bullet) in mass can be detected experimentally with a minimum of difficulty and absolutely no a priori information about the structure. Damage detection algorithms without a priori knowledge of any model and without first estimating a model are potentially useful for informing maintenance personnel of the existence of damage in a variety of civilian and military applications. These applications range from civil structures and infrastructures to manufacturing product imperfections to helicopter components. The simplicity of the algorithm allows on board diagnostics when coupled with self-sensing piezoceramic embedded composites. ACKNOWLEDGEMENTS
The authors gratefully acknowledge the support of ARO grant number DAAL 03-92-G-0180 under the supervision of Dr Gary Anderson. REFERENCES 1. P. F. R and N. A 1990 Journal of Sound and Vibration 138, 381–388. Identification of crack location and magnitude in a cantilevered beam from the vibration modes. 2. C. D, P. Q. Z, W. Q. F and T. C. H 1994 Proceedings of the 12th IMAC, Honolulu, HW, 98–105. The sensitivity study of the modal parameters of a cracked beam. 3. H. T. B, D. J. I, D. J. L and Y. W 1996 Journal of Sound and Vibration 191, 859–880. An experimentally validated damage detection theory in smart structures. 4. J. C, G. R. T and K. W 1994 Proceedings of the 12th IMAC, Honolulu, HW, 778–786. A simplified approach to the numerical and experimental modelling of the dynamics of a cracked beam. 5. Y. Y, D. S. B and R. E. S 1993 ACC Proceedings, San Francisco, CA, 3021–3023. Frequency domain identification for robust large space structures control design. 6. D. J. I 1994 Engineering Vibration. Englewood Cliffs, NJ: Simon and Schuster. 7. D. K. L and G. K 1994 AIAA Proceedings, Hilton Head, SC, 192–198. Location and estimation of damage in a beam using different identification algorithms.
APPENDIX A
A first step in assessing the technique’s feasibility is to set up a simulation for the longitudinal vibration of a cantilevered bar. The bar is divided into ne elements of length si allowing each node one single degree of freedom [6]. The stiffness matrix for the ith element is Ki =
$
Ei Ai 1 si −1
%
−1 EA = i i K i 1 si
(A1)
where Ei is Young’s modulus and Ai the cross-sectional area. The mass matrix is Mi =
$ %
rAi si 2 6 1
1 rA s = i i M i 2 6
(A2)
where r is the mass density. Since each element is attached at the ith and (i + 1)th nodes,
. . .
420
the global matrices can be obtained by direct assembly. Considering the algorithms provided in [7] the expanded stiffness matrix in global co-ordinates is Khi =
Ei Ai Bi K i BiT si
(A3)
introducing the matrix Bi Bi = [ei ei + 1 ]
(A4)
where ej is the jth unit vector, i.e. all elements are zero except for a value 1 in the jth position. Summing the individual stiffness matrices results in ne
ne Ei Ai Bi K i BiT = s Khi . si i=1 i=1
Kh = s
(A5)
The global mass matrix is assembled in exactly the same manner with rAi si Bi M i BiT 6
(A6)
ne rAi si Bi M i BiT = s Mhi . 6 i=1 i=1
(A7)
Mhi = which yields ne
Mh = s
The subscript h denotes quantities from the healthy model. The damage is represented by changes of the cross-sectional area which equal losses of stiffness and mass. While this is not the most advanced model of damage, it does provide a reasonable set of data in a controlled environment for use in testing the approach. Note that such a model is just a subset of all damage models. The cross-section for the ith element Ai will be varied by the scalar factor ai (1 − ai )Ai for 0 E ai E 1
(A8)
which is called the damage coefficient of the ith element of the beam. Suppose that any element may have sustained damage as sketched in Fig. A1. The damaged global stiffness and mass matrices then will change to ne
ne Ei (1 − ai )Ai Bi K i BiT = s (1 − ai )Kdi si i=1 i=1
Kd = s
Figure A1. Bar with reduction in cross-sectional area.
(A9)
ne ne r(1 − ai )Ai si Md = s Bi M i BiT = s (1 − ai )Mdi 6 i=1 i=1
421 (A10)
The subscript d denotes quantities from the damaged model. The beam model was set up to develop further the results obtained by axial excitations. Referring to Euler–Bernoulli beam elements, stiffness and mass matrices are given by equations (A11) and (A12). All constants are defined as before and Ii denotes the second moment of inertia. 6si −12 6si L K 12 2 G 6si 4si −6si 2si2 G Ei Ii Ki = 3 G G = Esi3Ii K i si −12 −6si 12 −6si i G G k 6si 2si2 −6si 4si2 l
(A11)
22si 54 −13si L K 156 G 22si 4si2 13si −3si2 G rA s i si Mi = i iG G = rA M . 420 54 13si 156 −22si 420 i G G k −13si −3si2 −22si 4si2 l
(A12)
Just as in Section 2.1, the global stiffness and global mass matrices are given by equations (A5) and (A7), respectively. The transformation matrix appears a little different from (A4) due to two degrees of freedom at each node. In particular Bi = [e2i − 1 e2i
e2i + 1
e2i + 2 ].
(A13)
The damage in a beam element is defined as a change in the flexural rigidity Ei Ii that will not affect the mass of the beam. It is assumed also that only sustained damage is modelled. Following the example of the bar, a damage coefficient decreases the values of the corresponding matrix elements. In this case, only the stiffness matrix is altered as indicated (1 − di )
Ei Ii K for 0 E di E 1. si3 i
(A14)
The global stiffness matrix then becomes ne
ne (1 − di )Ei Ii Bi K i BiT = s (1 − di )Kdi 3 si i=1 i=1
Kd = s
(A15)
and the global mass matrix remains ne ne rAi si Bi M i BiT = s Mdi = s Mhi . 420 i=1 i=1 i=1 ne
Md = Mh = s
(A16)
Suppose the damage is modeled as a crack such as a cut with depth Dhi shown in Fig. A2. In that case, it is obvious that of all constants involved in equation (A11) the
. . .
422
Figure A2. Bar with reduction in cross-sectional area.
damage coefficient di to represent the crack it is possible to derive a correlation between the depth of the crack and the percentage of stiffness reduction. Thus Ixx (Dhi ) = (1 − di )Ixx .
(A17)
The second moment of the area of the cracked element about the x-axis is denoted as Ixx (Dhi ) =
g
h/2 − Dhi
y 2bdy
−h/2
= 121 bh 3 − 14 bh 2Dhi + 12 bhDhi2 − 13 bDhi3
(A18)
and of the area of the healthy elements Ixx = 121 bh 3,
(A19)
respectively. Using equations (A18) and (A19) in (A17) results in
0 1 0 1 0 1 Dhi h
3
− 32
Dhi h
2
+ 34
Dhi − 14 di = 0. h
(A20)
Since the ratio Dhi /h is proportional to the loss of mass of the cracked element, equation (A20) evaluates the correlation of stiffness reduction and mass reduction for this one particular element. In Fig. A3 the relation between the two damage ratios di and Dhi /h is plotted. The assumption that damage, as modeled before, will not alter the mass but stiffness exclusively, needs to be applied carefully. It is shown that there is an
Figure A3. Damage ratio correlation.
423
immediate interrelation between the two quantities when damage is modeled as suggested in Fig. A2. Due to this result the global mass matrix defined in equation (A16) will now become ne
Md = s i=1
ne
=s i=1
0
1−
Dhi rAi si B M B T h 420 i i i
1
0
1−
Dhi Mdi with Dhi = f(di ) h
1
(A21)
APPENDIX B
The dynamic model for both undamaged mechanical structures, beam and bar is given by Mh u¨ (t) + Cu˙ (t) + Kh u(t) = F(t).
(B1)
The damping given in equation (B1) is assumed to be proportional, depending on both global stiffness and global mass matrices of the form C = aMh + bKh .
(B2)
Here, the constants a and b may be defined freely, yet it remains difficult to understand how the proportional damping will actually affect the model’s response. The damping ratio zj =
a bv + j 2vj 2
(B3)
takes remedial action. It determines the rate of decay of the jth natural frequency. Using expression (B3) derived in [6] in equation (B2) leads to a = 2v1 v2 b=2
z1 v2 − z2 v1 v22 − v12
z2 v2 − z1 v1 v22 − v12
(B4) (B5)
where the first and second natural frequencies are chosen, in regard to their significant influence on the system’s response. Substituting equation (B4) and (B5) into (B2) provides a more transparent damping behaviour. In the case of damage the dynamic behaviour will be modeled by Md u¨ (t) + Cu˙ (t) + Kd u(t) = F(t).
(B6)
Damping will be calculated as for the undamaged system but using Kd and Md , the global matrices of the damaged model. Note again that these models are needed to provide data and are not part of the ‘detection’ algorithm proposed here.