- Email: [email protected]
Numerical investigation of a new damage detection method based on proper orthogonal decomposition
ARTICLE IN PRESS Mechanical Systems and Signal Processing Mechanical Systems and Signal Processing 21 (2007) 1346–1361 www.elsevier.com/locate/jnlabr/ymssp
Numerical investigation of a new damage detection method based on proper orthogonal decomposition Ugo Galvanetto, George Violaris1 Department of Aeronautics, Imperial College London, Prince Consort Road, London SW7 2BY, UK Received 8 July 2005; received in revised form 9 December 2005; accepted 14 December 2005 Available online 8 February 2006
Abstract The paper presents a new idea for detecting and locating damage in structures by applying concepts derived from the theory of proper orthogonal decomposition. The idea is investigated by simulating tests on two beams and provides promising results. The effect of noise on the performance of the detection method is taken into account in the simulations. The paper paves the way to an experimental application of the newly proposed method. r 2006 Elsevier Ltd. All rights reserved. Keywords: Damage detection; Damage location; Proper orthogonal decomposition
1. Introduction Engineering systems in different fields, civil, mechanical, aeronautical, etc. are in continuous use despite ageing and the associate possibility of damage growth. Damage is usually defined in terms of a comparison between two different states of the system, one of which is generally assumed to represent the initial, undamaged state. The ability of monitoring the structural health of structural systems is becoming increasingly important. Non-destructive testing indicates the possibility of examining the integrity of a structure using a technique that does not risk introducing damage in the structure itself and is usually used for the detection of defects concentrated in small portions of the structure. A complete review of the different NDT methods is given in [1]. Many NDT methods are highly effective when applied locally, i.e. to a small portion of the structure, so that the examination of the whole structure requires several applications of the non-destructive technique. It would be useful to have global NDT methods that could, quickly and cheaply, reveal the presence of damage in an area of the structure to which a local technique could then be effectively applied for a more precise assessment. Damage identification based on variations of vibration properties is one of the few methods that monitor the changes in the whole structure simultaneously [2,3]. Although the basis of vibration based damage detection appears intuitive, its actual application poses many significant challenges. The most fundamental challenge is Corresponding author. Tel.: +44 (0) 20 7594 5150. 1
E-mail address: [email protected] (U. Galvanetto). Post-graduate student.
0888-3270/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2005.12.007
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1347
the fact that damage is typically a local phenomenon and may not significantly influence the low-frequency global response of a structure, which is usually measured during vibration tests. Vibration characteristics are those determined from the vibration responses of a structure. These are mainly the natural frequencies and the corresponding mode shapes and damping factors. The natural modes are the vibration property which has to be monitored in order not only to detect but also to locate the presence of damage. A survey on the use of natural frequency changes for damage detection is presented in [4] where it is concluded that the shift in natural frequencies has some important practical limitations. Modal shape changes seem to be more effective: in [5] the Laplacian operator on the mode shapes is used and in [6] changes in the mode shape curvature are monitored to detect damage. Both methods seem to work properly with severe damage but are too sensitive to the presence of noise. Despite the fact that several of the proposed vibration methods perform well in numerical simulations, their performances vary when put on trial in experimental test cases in such a way that no definite conclusion on their practical reliability can be drawn [7]. This paper presents the idea of characterising the vibration properties of a structure, and therefore their changes, by means of the proper orthogonal decomposition (POD). The POD, also known as Karhunen-Loe`ve decomposition [8], is emerging as a powerful experimental tool in dynamics and vibrations. It ‘provides a basis for the modal decomposition of an ensemble of functions, such as data obtained in the course of experiments. Its properties suggest that it is the preferred basis to use in various applications. The most striking of these is optimality: POD provides the most efficient way of capturing the dominant components of an infinite-dimensional process with only finitely many, and often surprisingly few, modes’ [8]. The Proper Orthogonal Modes (POMs) capture more energy per identified mode than any other orthogonal complete basis and can be identified with ease for non-linear and non-smooth systems. Moreover, the energy distribution between POMs, which is defined by the corresponding POV, can help identifying the most important modes [8–10]. The method may be used to quantify spatial coherence in turbulence problems [8–10] and oscillating structures monitored by several sensors [11] and determine the number of active state variables in a dynamic system. Extensive applications of POD to structural dynamics were carried out by Feeny et al. [12,13], Vakakis and co-workers [14–16] and Golinval and co-workers [17]. In the present work the POD theory will be applied to assess the structural integrity of beams. All data processed with methods based on the POD have been generated numerically with a finite element code and, therefore this paper constitutes a first step towards the evaluation of the proposed damage detection method. Applications of the same method to data generated in real experiments will constitute a second step of this research. POD is closely related to the Singular Value Decomposition [18] and the Principal Component Analysis [19], methods which have already been used in damage detection. In the present work the POD methods are applied with a dense grid of sensors which measure displacement-like quantities and can therefore provide the characteristic shape of the deformed structure during the steady-state dynamics. The changes in such a shape are used to detect and locate the damage. Such approach is different from those adopted in [18] and [19]. The paper is organised as follows: Section 2 summarises the main ideas of the POD from an applied point of view. Section 3 describes the numerical model, the computational techniques and the numerical noise generator used to simulate the experimental set-up. Section 4 describes the choice of the damage indicator, Section 5 presents the results obtained with clean and noisy data. Finally Section 6 draws some preliminary conclusions. 2. Fundamental ideas on the POD The main objective of this paper is to present an idea to detect and locate damage in engineering structures by computing their POMs. Application of POD to vibrating structures typically requires the experimental acquisition of displacements (or accelerations) at N locations of a dynamic system. The recorded values of the displacements (or accelerations) at a given time t are labelled d 1 ðtÞ; d 2 ðtÞ; . . . d N ðtÞ. If the displacements are sampled for M times, one can form displacement history arrays, such that d i ¼ ðd i ðt1 Þ; d i ðt2 Þ; . . . d i ðtM ÞÞT for i ¼ 1; 2; . . . N. In performing the POD these displacement histories are normalised by subtracting the mean value d¯ i . The vectors ai are then formed as ai ¼ d i d¯ i 1,
(1)
ARTICLE IN PRESS 1348
U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
where 1 is a vector of dimension M whose components are all equal to unity. Vectors ai are then used to form an M N matrix: A ¼ ½a1 ; a2 ; . . . ; aN .
(2)
Each row in A represents the different positions of the recorded points of the structure at a particular instant in time, for that reason it is often called a snapshot. Each column in A, as previously said, represents the zeroaveraged time series of a single point of the structure. The correlation matrix is then formed: R ¼ ð1=MÞAT A:
(3)
R is real and symmetric, therefore its eigenvectors form an orthogonal basis. The eigenvectors of R are the POMs and the eigenvalues the POVs of the system. When applied to linear systems the optimality features of POD greatly reduce the amount of data to be considered. In fact in [17] it has been shown that: ‘the forced harmonic response of a linear system is captured by a single POM. Nevertheless all the linear natural modes (LNMs) are in general necessary to reconstruct the response. This property is independent on the mass distribution and underlines the optimality of the POMsy The convergence of the POM to a natural mode is not guaranteed. The POM appears as a combination of all natural modes’. Therefore the application of POD techniques to a linear system concentrates in a unique POM all the information which is contained in an infinite number of LNMs. POD can result very beneficial also when applied to non-linear cases, but, in general, in this case more than one POM will have to be considered. In the general non-linear case not all the POMs are significant since some of them could carry a negligible amount of energy. There is the need for a criterion to identify the dominant mode shapes from the rest. This is achieved by using the POVs, which provide a measure of the strength of participation (or energy) of the corresponding mode to the signal. If li indicates any POV and POVs are in decreasing order so that l1 4l2 4 4lN , the number of dominant POMs is usually selected so that the eigenvalues satisfy the following condition for the smallest integer p [12]: Pp li X0:9999. (4) Pi¼1 N i¼1 li Eq. (4) suggests that the changes in the distribution of the energy of the signal between the POVs could probably be a useful tool for identifying structural changes in the system. In the present paper this idea will be briefly discussed in Section 4 but will not be pursued further. The paper will investigate the changes affecting the POMs. The procedure presented in this paper, based on proper orthogonal modes, retains the main advantage of vibration techniques, i.e. monitors the changes in the structure on a global basis [2], moreover it constitutes a significant step forward with respect to more traditional vibration techniques [7], based on the linear natural modes, for two main reasons:
The ‘optimality’ of the POMs only requires the use of a few (only one in the case of linear systems) POMs whereas applications using the LNMs often require the evaluation of several modes. Damage, particularly in its initial phases, is a local phenomenon which does not significantly affect the low frequency natural modes of a vibrating structure, whereas its effects are more noticeable on the higher modes [7]. For that reason methods based on the use of LNMs face the problem of deciding how many natural modes have to be examined and which ones are the most relevant to detect the possible presence of damage. On the contrary the dominant POM is associated with the largest eigenvalue and therefore uniquely determined. Moreover such optimality seems to make the POMs more robust with respect to the unavoidable presence of noise [16]. The computation of POMs and POVs, based on Eqs. (1)–(4) and on an eigenproblem solution, is extremely simple and fast. Moreover it does not rely on linearity assumptions and therefore the proposed method can be applied to both linear and non-linear systems.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1349
3. Fem models and numerical experiments The fem models investigated in this paper are shown in Fig. 1: the vertical beam is hinged at the top end where also a stiff rotational spring is applied; the horizontal beam is simply supported at both ends. In both cases the beam is homogeneous and uniform, is defined by the properties of Table 1 and is discretised with 40 equal elements. Both systems are damped with a Rayleigh type damping, which is proportional to the initial _ where C is the structural damping matrix, x_ the vector of nodal velocities, stiffness matrix: Cx_ ¼ ½aM þ bK0 x, M the structural mass matrix, K0 the structural initial stiffness matrix, a ¼ 0, b ¼ 4.38699 107 are the coefficients of Rayleigh damping. Large strains and shear deformation are taken into account by a Reissner type element [20]. A sinusoidal motion is imposed to the support of the cantilever beam (CB) and to either the supports or another node of the simply supported beam (SSB), in each case in the direction perpendicular to the axis of the beam. The dynamics of the beams is numerically integrated with a ‘mid-point algorithm’, a particular form of finite difference, which is very stable for long-term integrations [21]. For both models the POVs and POMs are computed for two cases: the undamaged beam and a damaged beam where one element has the Young’s modulus and shear modulus reduced with respect to the original value. In the case of the vertical CB the variables used to compute the POMs are the horizontal components of the nodal displacements. In the case of the horizontal SSB the variables used to compute the POMs are the vertical components of the nodal displacements. When the motion is applied to the supports the displacement components used to compute the POMs are relative to the points to which the motion is imposed. Only steady-state data are used in the computation, the initial transient is discharged. The results presented in the
(a)
(b) Fig. 1. Models of the beams: (a) cantilever beam with earth spring and its fem mesh, and (b) simply supported beam and its fem mesh. Table 1 Beam data Length: Width, b Thickness, h Cross-sec. area, A ¼ bh Second moment of area (bh3 =12): Density, r Young’s modulus, E Shear modulus, G
0.1880 m 0.00950 m 0.00023 m 2.185 106 m2 9.65 1015 m4 7840 kg/m3 2.0 105 MPa 1.4286 105 MPa
ARTICLE IN PRESS 1350
U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
next section were obtained by recording the nodal displacement values 8948 times. The time interval between two successive recorded data (sampling time Dt) was much bigger than the time step of the numerical integration but small enough to allow several recorded data per period of the imposed motion. Moreover the ratio (period of imposed motion/sampling time) was not a rational number so that a shift in the phase angle characterises data recorded in successive periods of the forcing motion. Using the terminology introduced above, the parameters have the following values: N ¼ 40 and M ¼ 8948. The forty vectors ai have 8948 components and R is a (40 40) square matrix with forty real eigenvalues and eigenvectors. The vectors ai are obtained from a FEM code, then with simple matrix operations the matrix R is computed, and with a mathematical program like Matlab, the eigenvalues and eigenvectors of R are easily obtained. The presence of damage and its location will be revealed by comparing the POMs of the original structure with those of the damaged structure. Extensive numerical investigations have shown that the numerical values of POM components and POVs are robust with respect to the two numerical parameters involved in their computation: sampling time (Dt) and number of time samples (M). If Dt and M are chosen in such a way that the motion of the system is described in sufficient detail no appreciable variation in the values of POMs and POVs can be noticed even for moderate variations of Dt and M. A few considerations about the usual linear vibration properties of the beams will be useful to locate the chosen frequency of the imposed motion among the natural frequencies of the investigated mechanical systems. The usual linear Bernoulli beam theory applied to a cantilever perfectly clamped beam with the mechanical data of Table 1 provides the following linear natural frequencies [22]: f1 ¼ 5.3, f2 ¼ 33.3, f3 ¼ 93.2 Hz. In a similar way linear natural frequencies for a SSB are given next: f1 ¼ 14.9, f2 ¼ 59.6, f3 ¼ 134.1 Hz. This information will be used in the next section to justify the choice of the frequency of the forcing. The size of the beam and the numerical values of the material constants do not have any particular relevance. The mathematical models of Fig. 1 were chosen only because already available, since they had been used in a previous experimental and numerical study [23]. Also the value of the damping parameter b was chosen because it provides a realistic behaviour. Other simulations, described in [24,25], carried out with different FEM programs and different structural models generated results analogous to those presented in Section 5. One of the main drawbacks of vibration methods, when applied to real situations, is their sensitivity to noise: often the changes in the vibration properties due to the damage in the structures are hidden by the high level of noise. It was decided to add white Gaussian noise to the time series generated by the finite element code. The Matlab command awgn was used to add white Gaussian noise to the displacement vectors di before they were processed to generate the correlation matrix R. The command awgn requires the specification of a scalar quantity snr (signal-to-noise ratio) that defines the ‘amplitude’ of the noise with respect to that of the clean signal. When in the next section the noise level is given by a particular value of snr it means that a noisy signal with such an snr has been added to the time series of each node. Moreover the noisy sequences affecting different nodes are uncorrelated, in this way severe experimental conditions were simulated. Fig. 2 gives a graphical description of the meaning of the noise level: the same signal, represented by the black line, is affected by three different levels of noise characterised by the values snr ¼ 20, 10, 1. The noisy signal is given by the grey line. 4. Choice of the damage indicator The main idea of the present work consists in choosing the changes in POMs as damage indicator. The POMs are not only a property of the structure and its boundary conditions, as the LNMs, but they depend also on the dynamics of the system. If a sinusoidal motion is applied to a structure its POMs change with the frequency of the excitation. It is expected that they should be similar to the LNMs when the frequency of the imposed motion is close to a natural frequency of the structure. For the CB simulations were carried out for 5.3 Hz (close to Mode I of the clamped CB), 8.5, 18.5 and 20.5 Hz (as steps between Mode I and II behaviour), 33.3 Hz (Mode II), 60.5 Hz (as a step between Mode II and III behaviour) and 90.5 Hz (close to Mode III). The POM shapes of the dominant POM are shown in Fig. 3. The figure confirms that the dominant POM visually
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1351
Matlab snr 20
1.5 1 0.5 0 1
13 25 37 49 51 73 85 97 109 121 133 145 157 169 181
-0.5 -1 -1.5 Signal
Noisy Signal
Matlab snr 10 2 1.5 1 0.5 0
1
13 25 37 49 51 73 85 97 109 121 133 145 157 169 181
-0.5 -1 -1.5 Signal
Noisy Signal
Matlab snr 1 3 2 1 0
1
13 25 37 49 51 73 85 97 109 121 133 145 157 169 181
-1 -2 -3 Signal
Noisy Signal
Fig. 2. Examples of different levels of noise.
resembles the LNM corresponding to a linear natural frequency close to that of the imposed motion. Moreover the first POV is several orders of magnitude larger than the others which means that the system under investigation can be approximated as a linear system even if the numerical program describes geometrically non-linear effects. In this respect the simulated results are similar to an experimental case where non-linearities may be present but are not activated by the mechanics of the system. Similar FE simulations were repeated for the SSB with motion imposed to the supports. Different frequencies were applied: 8.5 Hz (below Mode I), 20.5 Hz (above Mode I), 60.5 Hz (Mode II), 95.0 Hz (as a step between Mode II and III) and 134.5 Hz (Mode III). The Mode Shapes are shown in Fig. 4, which confirms the findings previously described for the CB.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1352 0.4
5.3Hz
8.5Hz
18.5Hz
20.5Hz
33.3Hz
60.5Hz
90.5Hz
Undamaged Beam Eigen Vectors
0.3
0.2
0.1
0 0
5
10
15
20
25
30
35
40
-0.1
-0.2
-0.3
Element
Fig. 3. Dominant POM for the cantilever beam (CB) at different frequencies of the imposed motion.
0.4 8.5Hz
20.5Hz
60.5Hz
95.0Hz
134.5Hz
Undamaged Beam Eigen Vectors
0.3
0.2
0.1
0 0
5
10
15
20
25
30
35
40
-0.1
-0.2
-0.3
Element
Fig. 4. Dominant POM for the simply supported beam (SSB) at different frequencies of the imposed motion.
The amplitude of the applied sinusoidal motion does not seem to have a considerable impact on the POMs of the beams. Several motions with different small amplitudes were applied and no appreciable difference in the POMs was noticed. All results presented in this paper were obtained with a displacement amplitude of 3 mm, which is the largest amplitude used in our computations.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1353
The POMs of a beam in its undamaged and damaged states are computed and presence and location of damage are inferred from the changes induced in the POMs by the damage. The principal hypothesis on which the method is based is that the system is linear in both states, undamaged and damaged. Such an assumption is true in the simulations, where damage is represented by a change in the material constants of one element, but the linear elastic nature of the constitutive law is not affected and non-linear geometrical effects are taken into consideration in the FEM algorithm but their magnitude is negligible. In real damaged structures the linearity assumption is probably verified in the case of initial small damage or also for more severe damage if the amplitude of the imposed oscillations is small. The important consequence of the assumed linearity is that the steady-state dynamics under harmonic excitation is captured by a unique POM, in accordance with the previously mentioned theory of [17]. Therefore the dynamic behaviour of the system is fully captured by only one POM which contains in its variations due to damage the information about damage that will be used in the present work. Moreover since the POM is defined along the axis of the beam its changes may also indicate the location of the damaged point of the structure. In the numerical models damage is introduced in an element (often element 20) as a reduction in Young’s and shear moduli. Then the same motions used to compute the dominant POMs in the undamaged structures, those shown in Figs. 3 and 4, are imposed to the damaged beams. In the next paragraph several figures will be plotted to show how the method works. The quantity plotted along the vertical axes of Figs. 5 and 6 (and of many other figures in the following pages) is the difference DPOM defined as DPOMðiÞ ¼ POMu ðiÞ POMd ðiÞ;
i ¼ 1; 2; . . . ; 40,
(5)
where POMu(i) is the value of the dominant POM in the undamaged case at the ith node of the mesh and POMd(i) is the value of the dominant POM in the damaged case at the same node. The horizontal axis of the figures represents the length along the axis of the beam measured in number of elements of the fem mesh. As previously noticed the changes in the POVs could be another candidate as damage indicator, but the variations of the POVs due to damage do not seem easy to interpret, moreover they do not provide any information about the possible location of damage. For these reasons DPOM rather than DPOV has been chosen as damage indicator.
0.0008
Delta Eigen Vector
-0.0002 0
5
10
15
20
25
30
33.3Hz
60.5Hz
35
-0.0012
-0.0022
-0.0032
-0.0042
Element 5.3Hz
8.5Hz
18.5Hz
20.5Hz
Fig. 5. DPOM at different frequencies of the imposed motion for the CB.
90.5Hz
40
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1354 0.0078
8.5Hz
20.5Hz
60.5Hz
95.0Hz
134.5Hz
Delta Eigen Vector
0.0058
0.0038
0.0018
-0.0002
0
5
10
15
20
25
30
35
40
-0.0022
-0.0042
Element
Fig. 6. DPOM at different frequencies of the imposed motion for the SSB.
5. Results 5.1. Damage detection and location The differences between the POMs of the original structures and those of the damaged ones are shown in Figs. 5 and 6. In both figures the plot of the difference has a sudden change in slope at the point where damage is located, element 20. Damage seems to be clearly detectable at any frequency. We observe that the simple examples shown in the present work illustrate some potentialities of the new method but are clearly not sufficient to demonstrate its effectiveness, especially for more complex structures examined in situ. The fact that the damage seems to be detectable at any frequency might be due to the simplicity of the examined cases. Usually damage in the early stages affects the higher LNMs and if the excitation frequency is far from those of such higher modes, little damage information may be measured and the damage may become difficult to detect. All computations presented in the remainder of the paper were obtained with a frequency of 8.5 Hz, which is above f1 for the CB and below f1 for the SSB. The main steps of the method proposed in the paper are given in the block diagram of Fig. 7.
5.2. Damage severity The mechanical properties of element 20 were reduced by 10%, 25%, 40% and 50% (see Table 2) to simulate increasing damage severity. The diagrams of the corresponding DPOM are plotted in Figs. 8 and 9 which show that more severe damage generates bigger slope discontinuities in the DPOMs. Exciting the SSB by moving its two supports would not be practical in an experimental environment. For that reason in the next examples the motion will be applied to a node of the beam model. It was found that exciting close to the centre of the beam produced results most similar to exciting through the supports and this is probably due to the almost symmetry of the beam. However, due to the asymmetric support conditions the behaviour of the model is not fully symmetric and it was found that exciting closer to node 0 (than to node 40) would produce the
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1355
Table 2 Mechanical Properties of the damaged beam
Undamaged 10% Damage 25% Damage 40% Damage 50% Damage
Young’s modulus (GPa)
Shear modulus (GPa)
Density (g/cm3)
200 180 150 120 100
143 129 108 86 72
7.84 7.84 7.84 7.84 7.84
1) Apply excitation to undamaged structure.
2) Compute dominant POM of the undamaged structure.
3) Apply excitation to damaged structure.
4) Compute dominant POM of the damaged structure.
5) Compute the difference of the POMs. Damage will be located where the slope of the difference has an evident discontinuity.
Fig. 7. Block diagram describing the steps of the POD damage detection method.
clearest results for detecting the presence of damage in any node. The plots shown in the next sections were obtained by applying the external motion to node 11 of the SSB.
5.3. Different damage positions The effectiveness of the method at locating damage can vary with the position of the damage in the structure. Damage can be introduced in different elements. The effect of this is shown in Figs. 10 and 11, which were generated by introducing a 25% reduction in the mechanical properties in different elements of the mesh. Fig. 10 suggests that if damage is in certain zones of the structure it could be difficult to detect. In particular when damage is located at node 35 of the CB (or closer to the free end) the diagram does not seem to indicate the presence of damage in a clear way. The method seems to work well irrespective of where the damage is located in the SSB.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1356
0.0014
Delta Eigen Vector
0.0009
0.0004
-0.0001 0
5
10
15
20
25
30
35
40
-0.0006
Damage 10% Damage 25% Damage 40% Damage 50%
-0.0011
-0.0016
Element
Fig. 8. Effect of damage severity, CB, damage in element 20.
0.003
Damage 10% 0.0025
Damage 25% Damage 50%
0.002
Delta Eigen Vector
0.0015 0.001 0.0005 0 0
5
10
15
20
25
30
35
40
-0.0005 -0.001 -0.0015 -0.002
Element
Fig. 9. Effect of damage severity, SSB, damage in element 20. Motion imposed to nodes 1 and 41.
5.4. Noise White Gaussian noise was added to previous results of the CB and of the SSB. Figs. 12 and 13 illustrate that as the severity of damage is increased, the tolerance of the newly proposed method to noise also increases. In
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1357
0.00038
-0.00002
0
5
10
15
20
25
30
35
40
-0.00022
-0.00042
-0.00062
Element Damage Node 5 Damage Node 25
Damage Node 10 Damage Node 30
Damage Node 15 Damage Node 35
Damage Node 20
Fig. 10. Twenty-five per cent damage at various locations of the CB.
Damage Node 11
Damage Node 20
Damage Node 31
0.00125
0.00075
Delta Eigen Vector
Delta Eigen Vector
0.00018
0.00025
0
5
10
15
20
25
30
35
-0.00025
-0.00075
Element
Fig. 11. Twenty-five per cent damage at various locations of the SSB, excitation at node 11.
40
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1358 0.0003
0.0002 Effect of noise (10% damage at element 20)
Delta Eigen Vector
0.0001
Noise level 1%
0 0
5
10
15
20
25
30
35
40 Noise level 5%
-0.0001
Noise level 10% Noise level 20%
-0.0002
-0.0003
-0.0004
Element
0.0008 Effect of noise (25% damage at element 20)
0.0006
Noise level 1%
Delta Eigen Vector
0.0004
Noise level 5% 0.0002 0
Noise level 10% 0
5
10
15
20
25
30
35
40
Noise level 20% Noise level 30%
-0.0002 -0.0004 -0.0006 -0.0008
Element
0.0015 Effect of noise (50% damage at element 20)
0.001
Noise level 1%
Delta Eigen Vector
0.0005
Noise level 5% 0
0
5
10
15
20
25
30
35
40
Noise level 10% Noise level 20%
-0.0005 Noise level 30% Noise level 40%
-0.001
Noise level 50%
-0.0015
-0.002 Element
Fig. 12. Effect of noise on damage detection in the CB.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1359
Table 3 Relationship between % value and snr value of the noise level % value
1%
5%
10%
20%
30%
40%
50%
snr
48.417
34.250
28.148
22.046
18.477
15.945
13.980
Figs. 12 and 13 the noise level is indicated by a percent value which is linked to the snr parameter of Matlab according to Table 3. For 10% damage the noise dominates if its percent value is above 10%, in both the CB and the SSB. For 25% damage the tolerance is increased and damage can be seen clearly even at 20% noise level in the SSB, whereas its detection is more problematic in the CB. However, the signal becomes noise dominated for higher levels of noise. Finally, for 50% damage, noise has limited effect on the process of damage detection in both models. While the noise will distort the Mode shapes its effect is small compared to the magnitude of the POM differences.
5.5. Remarks Possibly one of the main drawbacks of the POD damage detection method, as presented in this paper, is that it is based on data coming from a dense grid of sensors distributed on the structure. That would make it expensive and impractical, especially for large structures. The authors believe that the most likely applications of the method under investigation would probably be in the quality control process of industrial production and would allow for a fast and safe detection of structural components (blades, stiffened panels, tubes, etc.) which do not meet the required standards. In such cases it would probably be possible to compute the POMs using data obtained with non-contact optical devices. Another possible field of application could be represented by smart structures already equipped with embedded sensors and actuators which would be used to implement the new technique. Moreover the method could be used also for real-time damage detection, for example in dynamic fatigue tests with no interruption of the test itself [26].
6. Conclusions A new idea for detecting and locating damage in structures has been presented. It is based on concepts of the theory of POD and has been investigated at a numerical level for beams.
Damage is revealed by a sudden change in the slope of the curve showing the value of DPOM, the difference between the values of the dominant POM of the undamaged and damaged structures. An increase in the intensity of damage generates an increase in the value of the change in slope of DPOM. It can be problematic to locate the damage if it affects certain zones of the structure, but in the investigated examples that often happens in the region where stresses and strains are the smallest, and therefore where there is a smaller probability of having damage. The presence of noise can reduce the capability of the method to locate the damage.
The methodology seems promising since it is capable of locating the damage and presents a degree of robustness with respect to the presence of noise. The method is attractive because it is global, easy to implement and can be extended to structures of more complicated shape.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1360 0.0004 0.0003
Effect of noise (10% damage at element 20)
Delta Eigen Vector
0.0002 0.0001
Noise level 1% 0
0
5
10
15
20
25
30
35
40
Noise level 5% Noise level 10%
-0.0001
Noise level 20% -0.0002 -0.0003 -0.0004 Element 0.001
Effect of noise (25% damage at element 20)
Delta Eigen Vector
0.0005
Noise level 1% Noise level 5%
0 0
5
10
15
20
25
30
35
40
Noise level 10% Noise level 20%
-0.0005
-0.001 Element Effect of noise (50% damage at element 20)
0.0024 0.0019
Noise level 1%
Delta Eigen Vector
0.0014
Noise level 5% 0.0009 Noise level 10% 0.0004 -0.0001 0
Noise level 20% 5
10
15
20
25
30
35
40
Noise level 30% Noise level 40%
-0.0006
Noise level 50% -0.0011 -0.0016 -0.0021 Element
Fig. 13. Effect of noise on damage detection in the SSB.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1361
References [1] S. Ness, C.N. Sherlock (Eds.), Nondestructive testing overview, ASTN Nondestructive testing Handbook, Vol. 10, American Society of NDT, 1996. [2] S.W. Doebling, C.R. Farrar, M.B. Prime, A summary review of vibration-based damage identification methods, The Shock and Vibration Digest 30 (2) (1998) 91–105. [3] P. Cawley, R.D. Adams, The location of defects in structures from measurements of natural frequencies, Journal of Strain Analysis 14 (1979) 49–57. [4] O.S. Salawu, Detection of structural damage through changes in frequency: a review, Engineering Structures 19 (9) (1997) 718–723. [5] C.P. Ratcliffe, Damage detection using a modified Laplacian operator on mode shape data, Journal of Sound and Vibration 204 (3) (1997) 505–517. [6] K. Pandey, M. Biswas, M.M. Samman, Damage detection from changes in curvature mode shapes, Journal of Sound and Vibration 145 (1991) 321–332. [7] Y.K. Ho Structural damage identification using changes in vibration characteristics, Ph.D. thesis, Imperial College London, 2000. [8] P. Holmes, J.L. Lumley, G. Bekooz, Turbulence, coherent structures, dynamical systems and symmetry, Cambridge University Press, Cambridge, 1996. [9] J.L. Lumley, Stochastic Tools in Turbulence, Academic Press, New York, 1970. [10] L. Sirovich, Turbulence and dynamics of coherent structures. Part 1: Coherent structures, Quarterly of Applied Mathematics 45 (3) (1987) 561–571. [11] J.P. Cusumano, M.T. Sharkady, B.W. Kimble, Experimental measurements of dimensionality and spatial coherence in the dynamics of a flexible beam, Philosophical Transactions of the Royal Society of London Series A 347 (1994) 421–438. [12] R.V. Kappagantu, B.F. Feeny, Part 2: Proper orthogonal modal modeling of a frictionally excited beam, Nonlinear Dynamics 23 (2000) 1–11. [13] B.F. Feeny, On proper orthogonal co-ordinates as indicators of modal activity, Journal of Sound and Vibration 255 (2002) 805–817. [14] M.F. Azeez, A.F. Vakakis, Proper orthogonal decomposition (POD) of a class of vibroimpact oscillations, Journal of Sound and Vibration 240 (5) (2001) 859–889. [15] X. Ma, A.F. Vakakis, L.A. Bergman, Karhunen-Loe`ve modes of a truss: Transient response reconstruction and experimental verification, AIAA Journal 39 (4) (2001) 687–696. [16] I.T. Georgiou, I. Schwartz, ‘POD (K-L) analysis of the interaction free dynamics of a beam/pendulum configuration: an experimental study’, Proceedings of DETC’01, Pittsburgh, PA, September 9–12, 2001. [17] G. Kerschen, J.C. Golinval, Physical interpretation of the proper orthogonal modes using the singular value decomposition, Journal of Sound and Vibration 249 (2002) 849–865. [18] R. Ruotolo, C. Surace, Using SVD to detect damage in structures with different operational conditions, Journal of Sound and Vibration 226 (1999) 425–439. [19] H. Soon, J.A. Czarnecki, C.R. Farrar, Structural health monitoring using statistical process control, Journal of Structural Engineering (2000) 1356–1363. [20] E. Reissner, On one-dimensional finite strain beam theory: The plane problem, Journal of Applied Mathematics and Physics (ZAMP) 23 (1972) 795–804. [21] P.B. Bornemann, U. Galvanetto, Discrete dynamics of implicit time integration schemes for a dissipative system, Communications in Numerical Methods in Engineering 19 (2003) 737–748. [22] R.R. Craig Jr., Structural Dynamics, Wiley, New York, 1981. [23] P.B. Bornemann, Time integration algorithms for the steady states of dissipative non-linear dynamic systems, Ph.D. thesis, Imperial College London, 2003. [24] G. Violaris, Damage detection in beam structures using the proper orthogonal decomposition, Master Thesis, Department of Aeronautics, The composites Centre, Imperial College London, 2004. [25] D. Hunter, Damage detection in isotropic, cantilever plates using Proper Orthogonal Decomposition, Final Year Project Report, Department of Aeronautics, Imperial College London, 2005. [26] M. Wiercigroch, C.H. Foong, W.F. Deans, An experimental rig to investigate fatigue crack growth under dynamic loading, Meccanica 38 (2003) 19–31.
Numerical investigation of a new damage detection method based on proper orthogonal decomposition Ugo Galvanetto, George Violaris1 Department of Aeronautics, Imperial College London, Prince Consort Road, London SW7 2BY, UK Received 8 July 2005; received in revised form 9 December 2005; accepted 14 December 2005 Available online 8 February 2006
Abstract The paper presents a new idea for detecting and locating damage in structures by applying concepts derived from the theory of proper orthogonal decomposition. The idea is investigated by simulating tests on two beams and provides promising results. The effect of noise on the performance of the detection method is taken into account in the simulations. The paper paves the way to an experimental application of the newly proposed method. r 2006 Elsevier Ltd. All rights reserved. Keywords: Damage detection; Damage location; Proper orthogonal decomposition
1. Introduction Engineering systems in different fields, civil, mechanical, aeronautical, etc. are in continuous use despite ageing and the associate possibility of damage growth. Damage is usually defined in terms of a comparison between two different states of the system, one of which is generally assumed to represent the initial, undamaged state. The ability of monitoring the structural health of structural systems is becoming increasingly important. Non-destructive testing indicates the possibility of examining the integrity of a structure using a technique that does not risk introducing damage in the structure itself and is usually used for the detection of defects concentrated in small portions of the structure. A complete review of the different NDT methods is given in [1]. Many NDT methods are highly effective when applied locally, i.e. to a small portion of the structure, so that the examination of the whole structure requires several applications of the non-destructive technique. It would be useful to have global NDT methods that could, quickly and cheaply, reveal the presence of damage in an area of the structure to which a local technique could then be effectively applied for a more precise assessment. Damage identification based on variations of vibration properties is one of the few methods that monitor the changes in the whole structure simultaneously [2,3]. Although the basis of vibration based damage detection appears intuitive, its actual application poses many significant challenges. The most fundamental challenge is Corresponding author. Tel.: +44 (0) 20 7594 5150. 1
E-mail address: [email protected] (U. Galvanetto). Post-graduate student.
0888-3270/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2005.12.007
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1347
the fact that damage is typically a local phenomenon and may not significantly influence the low-frequency global response of a structure, which is usually measured during vibration tests. Vibration characteristics are those determined from the vibration responses of a structure. These are mainly the natural frequencies and the corresponding mode shapes and damping factors. The natural modes are the vibration property which has to be monitored in order not only to detect but also to locate the presence of damage. A survey on the use of natural frequency changes for damage detection is presented in [4] where it is concluded that the shift in natural frequencies has some important practical limitations. Modal shape changes seem to be more effective: in [5] the Laplacian operator on the mode shapes is used and in [6] changes in the mode shape curvature are monitored to detect damage. Both methods seem to work properly with severe damage but are too sensitive to the presence of noise. Despite the fact that several of the proposed vibration methods perform well in numerical simulations, their performances vary when put on trial in experimental test cases in such a way that no definite conclusion on their practical reliability can be drawn [7]. This paper presents the idea of characterising the vibration properties of a structure, and therefore their changes, by means of the proper orthogonal decomposition (POD). The POD, also known as Karhunen-Loe`ve decomposition [8], is emerging as a powerful experimental tool in dynamics and vibrations. It ‘provides a basis for the modal decomposition of an ensemble of functions, such as data obtained in the course of experiments. Its properties suggest that it is the preferred basis to use in various applications. The most striking of these is optimality: POD provides the most efficient way of capturing the dominant components of an infinite-dimensional process with only finitely many, and often surprisingly few, modes’ [8]. The Proper Orthogonal Modes (POMs) capture more energy per identified mode than any other orthogonal complete basis and can be identified with ease for non-linear and non-smooth systems. Moreover, the energy distribution between POMs, which is defined by the corresponding POV, can help identifying the most important modes [8–10]. The method may be used to quantify spatial coherence in turbulence problems [8–10] and oscillating structures monitored by several sensors [11] and determine the number of active state variables in a dynamic system. Extensive applications of POD to structural dynamics were carried out by Feeny et al. [12,13], Vakakis and co-workers [14–16] and Golinval and co-workers [17]. In the present work the POD theory will be applied to assess the structural integrity of beams. All data processed with methods based on the POD have been generated numerically with a finite element code and, therefore this paper constitutes a first step towards the evaluation of the proposed damage detection method. Applications of the same method to data generated in real experiments will constitute a second step of this research. POD is closely related to the Singular Value Decomposition [18] and the Principal Component Analysis [19], methods which have already been used in damage detection. In the present work the POD methods are applied with a dense grid of sensors which measure displacement-like quantities and can therefore provide the characteristic shape of the deformed structure during the steady-state dynamics. The changes in such a shape are used to detect and locate the damage. Such approach is different from those adopted in [18] and [19]. The paper is organised as follows: Section 2 summarises the main ideas of the POD from an applied point of view. Section 3 describes the numerical model, the computational techniques and the numerical noise generator used to simulate the experimental set-up. Section 4 describes the choice of the damage indicator, Section 5 presents the results obtained with clean and noisy data. Finally Section 6 draws some preliminary conclusions. 2. Fundamental ideas on the POD The main objective of this paper is to present an idea to detect and locate damage in engineering structures by computing their POMs. Application of POD to vibrating structures typically requires the experimental acquisition of displacements (or accelerations) at N locations of a dynamic system. The recorded values of the displacements (or accelerations) at a given time t are labelled d 1 ðtÞ; d 2 ðtÞ; . . . d N ðtÞ. If the displacements are sampled for M times, one can form displacement history arrays, such that d i ¼ ðd i ðt1 Þ; d i ðt2 Þ; . . . d i ðtM ÞÞT for i ¼ 1; 2; . . . N. In performing the POD these displacement histories are normalised by subtracting the mean value d¯ i . The vectors ai are then formed as ai ¼ d i d¯ i 1,
(1)
ARTICLE IN PRESS 1348
U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
where 1 is a vector of dimension M whose components are all equal to unity. Vectors ai are then used to form an M N matrix: A ¼ ½a1 ; a2 ; . . . ; aN .
(2)
Each row in A represents the different positions of the recorded points of the structure at a particular instant in time, for that reason it is often called a snapshot. Each column in A, as previously said, represents the zeroaveraged time series of a single point of the structure. The correlation matrix is then formed: R ¼ ð1=MÞAT A:
(3)
R is real and symmetric, therefore its eigenvectors form an orthogonal basis. The eigenvectors of R are the POMs and the eigenvalues the POVs of the system. When applied to linear systems the optimality features of POD greatly reduce the amount of data to be considered. In fact in [17] it has been shown that: ‘the forced harmonic response of a linear system is captured by a single POM. Nevertheless all the linear natural modes (LNMs) are in general necessary to reconstruct the response. This property is independent on the mass distribution and underlines the optimality of the POMsy The convergence of the POM to a natural mode is not guaranteed. The POM appears as a combination of all natural modes’. Therefore the application of POD techniques to a linear system concentrates in a unique POM all the information which is contained in an infinite number of LNMs. POD can result very beneficial also when applied to non-linear cases, but, in general, in this case more than one POM will have to be considered. In the general non-linear case not all the POMs are significant since some of them could carry a negligible amount of energy. There is the need for a criterion to identify the dominant mode shapes from the rest. This is achieved by using the POVs, which provide a measure of the strength of participation (or energy) of the corresponding mode to the signal. If li indicates any POV and POVs are in decreasing order so that l1 4l2 4 4lN , the number of dominant POMs is usually selected so that the eigenvalues satisfy the following condition for the smallest integer p [12]: Pp li X0:9999. (4) Pi¼1 N i¼1 li Eq. (4) suggests that the changes in the distribution of the energy of the signal between the POVs could probably be a useful tool for identifying structural changes in the system. In the present paper this idea will be briefly discussed in Section 4 but will not be pursued further. The paper will investigate the changes affecting the POMs. The procedure presented in this paper, based on proper orthogonal modes, retains the main advantage of vibration techniques, i.e. monitors the changes in the structure on a global basis [2], moreover it constitutes a significant step forward with respect to more traditional vibration techniques [7], based on the linear natural modes, for two main reasons:
The ‘optimality’ of the POMs only requires the use of a few (only one in the case of linear systems) POMs whereas applications using the LNMs often require the evaluation of several modes. Damage, particularly in its initial phases, is a local phenomenon which does not significantly affect the low frequency natural modes of a vibrating structure, whereas its effects are more noticeable on the higher modes [7]. For that reason methods based on the use of LNMs face the problem of deciding how many natural modes have to be examined and which ones are the most relevant to detect the possible presence of damage. On the contrary the dominant POM is associated with the largest eigenvalue and therefore uniquely determined. Moreover such optimality seems to make the POMs more robust with respect to the unavoidable presence of noise [16]. The computation of POMs and POVs, based on Eqs. (1)–(4) and on an eigenproblem solution, is extremely simple and fast. Moreover it does not rely on linearity assumptions and therefore the proposed method can be applied to both linear and non-linear systems.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1349
3. Fem models and numerical experiments The fem models investigated in this paper are shown in Fig. 1: the vertical beam is hinged at the top end where also a stiff rotational spring is applied; the horizontal beam is simply supported at both ends. In both cases the beam is homogeneous and uniform, is defined by the properties of Table 1 and is discretised with 40 equal elements. Both systems are damped with a Rayleigh type damping, which is proportional to the initial _ where C is the structural damping matrix, x_ the vector of nodal velocities, stiffness matrix: Cx_ ¼ ½aM þ bK0 x, M the structural mass matrix, K0 the structural initial stiffness matrix, a ¼ 0, b ¼ 4.38699 107 are the coefficients of Rayleigh damping. Large strains and shear deformation are taken into account by a Reissner type element [20]. A sinusoidal motion is imposed to the support of the cantilever beam (CB) and to either the supports or another node of the simply supported beam (SSB), in each case in the direction perpendicular to the axis of the beam. The dynamics of the beams is numerically integrated with a ‘mid-point algorithm’, a particular form of finite difference, which is very stable for long-term integrations [21]. For both models the POVs and POMs are computed for two cases: the undamaged beam and a damaged beam where one element has the Young’s modulus and shear modulus reduced with respect to the original value. In the case of the vertical CB the variables used to compute the POMs are the horizontal components of the nodal displacements. In the case of the horizontal SSB the variables used to compute the POMs are the vertical components of the nodal displacements. When the motion is applied to the supports the displacement components used to compute the POMs are relative to the points to which the motion is imposed. Only steady-state data are used in the computation, the initial transient is discharged. The results presented in the
(a)
(b) Fig. 1. Models of the beams: (a) cantilever beam with earth spring and its fem mesh, and (b) simply supported beam and its fem mesh. Table 1 Beam data Length: Width, b Thickness, h Cross-sec. area, A ¼ bh Second moment of area (bh3 =12): Density, r Young’s modulus, E Shear modulus, G
0.1880 m 0.00950 m 0.00023 m 2.185 106 m2 9.65 1015 m4 7840 kg/m3 2.0 105 MPa 1.4286 105 MPa
ARTICLE IN PRESS 1350
U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
next section were obtained by recording the nodal displacement values 8948 times. The time interval between two successive recorded data (sampling time Dt) was much bigger than the time step of the numerical integration but small enough to allow several recorded data per period of the imposed motion. Moreover the ratio (period of imposed motion/sampling time) was not a rational number so that a shift in the phase angle characterises data recorded in successive periods of the forcing motion. Using the terminology introduced above, the parameters have the following values: N ¼ 40 and M ¼ 8948. The forty vectors ai have 8948 components and R is a (40 40) square matrix with forty real eigenvalues and eigenvectors. The vectors ai are obtained from a FEM code, then with simple matrix operations the matrix R is computed, and with a mathematical program like Matlab, the eigenvalues and eigenvectors of R are easily obtained. The presence of damage and its location will be revealed by comparing the POMs of the original structure with those of the damaged structure. Extensive numerical investigations have shown that the numerical values of POM components and POVs are robust with respect to the two numerical parameters involved in their computation: sampling time (Dt) and number of time samples (M). If Dt and M are chosen in such a way that the motion of the system is described in sufficient detail no appreciable variation in the values of POMs and POVs can be noticed even for moderate variations of Dt and M. A few considerations about the usual linear vibration properties of the beams will be useful to locate the chosen frequency of the imposed motion among the natural frequencies of the investigated mechanical systems. The usual linear Bernoulli beam theory applied to a cantilever perfectly clamped beam with the mechanical data of Table 1 provides the following linear natural frequencies [22]: f1 ¼ 5.3, f2 ¼ 33.3, f3 ¼ 93.2 Hz. In a similar way linear natural frequencies for a SSB are given next: f1 ¼ 14.9, f2 ¼ 59.6, f3 ¼ 134.1 Hz. This information will be used in the next section to justify the choice of the frequency of the forcing. The size of the beam and the numerical values of the material constants do not have any particular relevance. The mathematical models of Fig. 1 were chosen only because already available, since they had been used in a previous experimental and numerical study [23]. Also the value of the damping parameter b was chosen because it provides a realistic behaviour. Other simulations, described in [24,25], carried out with different FEM programs and different structural models generated results analogous to those presented in Section 5. One of the main drawbacks of vibration methods, when applied to real situations, is their sensitivity to noise: often the changes in the vibration properties due to the damage in the structures are hidden by the high level of noise. It was decided to add white Gaussian noise to the time series generated by the finite element code. The Matlab command awgn was used to add white Gaussian noise to the displacement vectors di before they were processed to generate the correlation matrix R. The command awgn requires the specification of a scalar quantity snr (signal-to-noise ratio) that defines the ‘amplitude’ of the noise with respect to that of the clean signal. When in the next section the noise level is given by a particular value of snr it means that a noisy signal with such an snr has been added to the time series of each node. Moreover the noisy sequences affecting different nodes are uncorrelated, in this way severe experimental conditions were simulated. Fig. 2 gives a graphical description of the meaning of the noise level: the same signal, represented by the black line, is affected by three different levels of noise characterised by the values snr ¼ 20, 10, 1. The noisy signal is given by the grey line. 4. Choice of the damage indicator The main idea of the present work consists in choosing the changes in POMs as damage indicator. The POMs are not only a property of the structure and its boundary conditions, as the LNMs, but they depend also on the dynamics of the system. If a sinusoidal motion is applied to a structure its POMs change with the frequency of the excitation. It is expected that they should be similar to the LNMs when the frequency of the imposed motion is close to a natural frequency of the structure. For the CB simulations were carried out for 5.3 Hz (close to Mode I of the clamped CB), 8.5, 18.5 and 20.5 Hz (as steps between Mode I and II behaviour), 33.3 Hz (Mode II), 60.5 Hz (as a step between Mode II and III behaviour) and 90.5 Hz (close to Mode III). The POM shapes of the dominant POM are shown in Fig. 3. The figure confirms that the dominant POM visually
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1351
Matlab snr 20
1.5 1 0.5 0 1
13 25 37 49 51 73 85 97 109 121 133 145 157 169 181
-0.5 -1 -1.5 Signal
Noisy Signal
Matlab snr 10 2 1.5 1 0.5 0
1
13 25 37 49 51 73 85 97 109 121 133 145 157 169 181
-0.5 -1 -1.5 Signal
Noisy Signal
Matlab snr 1 3 2 1 0
1
13 25 37 49 51 73 85 97 109 121 133 145 157 169 181
-1 -2 -3 Signal
Noisy Signal
Fig. 2. Examples of different levels of noise.
resembles the LNM corresponding to a linear natural frequency close to that of the imposed motion. Moreover the first POV is several orders of magnitude larger than the others which means that the system under investigation can be approximated as a linear system even if the numerical program describes geometrically non-linear effects. In this respect the simulated results are similar to an experimental case where non-linearities may be present but are not activated by the mechanics of the system. Similar FE simulations were repeated for the SSB with motion imposed to the supports. Different frequencies were applied: 8.5 Hz (below Mode I), 20.5 Hz (above Mode I), 60.5 Hz (Mode II), 95.0 Hz (as a step between Mode II and III) and 134.5 Hz (Mode III). The Mode Shapes are shown in Fig. 4, which confirms the findings previously described for the CB.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1352 0.4
5.3Hz
8.5Hz
18.5Hz
20.5Hz
33.3Hz
60.5Hz
90.5Hz
Undamaged Beam Eigen Vectors
0.3
0.2
0.1
0 0
5
10
15
20
25
30
35
40
-0.1
-0.2
-0.3
Element
Fig. 3. Dominant POM for the cantilever beam (CB) at different frequencies of the imposed motion.
0.4 8.5Hz
20.5Hz
60.5Hz
95.0Hz
134.5Hz
Undamaged Beam Eigen Vectors
0.3
0.2
0.1
0 0
5
10
15
20
25
30
35
40
-0.1
-0.2
-0.3
Element
Fig. 4. Dominant POM for the simply supported beam (SSB) at different frequencies of the imposed motion.
The amplitude of the applied sinusoidal motion does not seem to have a considerable impact on the POMs of the beams. Several motions with different small amplitudes were applied and no appreciable difference in the POMs was noticed. All results presented in this paper were obtained with a displacement amplitude of 3 mm, which is the largest amplitude used in our computations.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1353
The POMs of a beam in its undamaged and damaged states are computed and presence and location of damage are inferred from the changes induced in the POMs by the damage. The principal hypothesis on which the method is based is that the system is linear in both states, undamaged and damaged. Such an assumption is true in the simulations, where damage is represented by a change in the material constants of one element, but the linear elastic nature of the constitutive law is not affected and non-linear geometrical effects are taken into consideration in the FEM algorithm but their magnitude is negligible. In real damaged structures the linearity assumption is probably verified in the case of initial small damage or also for more severe damage if the amplitude of the imposed oscillations is small. The important consequence of the assumed linearity is that the steady-state dynamics under harmonic excitation is captured by a unique POM, in accordance with the previously mentioned theory of [17]. Therefore the dynamic behaviour of the system is fully captured by only one POM which contains in its variations due to damage the information about damage that will be used in the present work. Moreover since the POM is defined along the axis of the beam its changes may also indicate the location of the damaged point of the structure. In the numerical models damage is introduced in an element (often element 20) as a reduction in Young’s and shear moduli. Then the same motions used to compute the dominant POMs in the undamaged structures, those shown in Figs. 3 and 4, are imposed to the damaged beams. In the next paragraph several figures will be plotted to show how the method works. The quantity plotted along the vertical axes of Figs. 5 and 6 (and of many other figures in the following pages) is the difference DPOM defined as DPOMðiÞ ¼ POMu ðiÞ POMd ðiÞ;
i ¼ 1; 2; . . . ; 40,
(5)
where POMu(i) is the value of the dominant POM in the undamaged case at the ith node of the mesh and POMd(i) is the value of the dominant POM in the damaged case at the same node. The horizontal axis of the figures represents the length along the axis of the beam measured in number of elements of the fem mesh. As previously noticed the changes in the POVs could be another candidate as damage indicator, but the variations of the POVs due to damage do not seem easy to interpret, moreover they do not provide any information about the possible location of damage. For these reasons DPOM rather than DPOV has been chosen as damage indicator.
0.0008
Delta Eigen Vector
-0.0002 0
5
10
15
20
25
30
33.3Hz
60.5Hz
35
-0.0012
-0.0022
-0.0032
-0.0042
Element 5.3Hz
8.5Hz
18.5Hz
20.5Hz
Fig. 5. DPOM at different frequencies of the imposed motion for the CB.
90.5Hz
40
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1354 0.0078
8.5Hz
20.5Hz
60.5Hz
95.0Hz
134.5Hz
Delta Eigen Vector
0.0058
0.0038
0.0018
-0.0002
0
5
10
15
20
25
30
35
40
-0.0022
-0.0042
Element
Fig. 6. DPOM at different frequencies of the imposed motion for the SSB.
5. Results 5.1. Damage detection and location The differences between the POMs of the original structures and those of the damaged ones are shown in Figs. 5 and 6. In both figures the plot of the difference has a sudden change in slope at the point where damage is located, element 20. Damage seems to be clearly detectable at any frequency. We observe that the simple examples shown in the present work illustrate some potentialities of the new method but are clearly not sufficient to demonstrate its effectiveness, especially for more complex structures examined in situ. The fact that the damage seems to be detectable at any frequency might be due to the simplicity of the examined cases. Usually damage in the early stages affects the higher LNMs and if the excitation frequency is far from those of such higher modes, little damage information may be measured and the damage may become difficult to detect. All computations presented in the remainder of the paper were obtained with a frequency of 8.5 Hz, which is above f1 for the CB and below f1 for the SSB. The main steps of the method proposed in the paper are given in the block diagram of Fig. 7.
5.2. Damage severity The mechanical properties of element 20 were reduced by 10%, 25%, 40% and 50% (see Table 2) to simulate increasing damage severity. The diagrams of the corresponding DPOM are plotted in Figs. 8 and 9 which show that more severe damage generates bigger slope discontinuities in the DPOMs. Exciting the SSB by moving its two supports would not be practical in an experimental environment. For that reason in the next examples the motion will be applied to a node of the beam model. It was found that exciting close to the centre of the beam produced results most similar to exciting through the supports and this is probably due to the almost symmetry of the beam. However, due to the asymmetric support conditions the behaviour of the model is not fully symmetric and it was found that exciting closer to node 0 (than to node 40) would produce the
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1355
Table 2 Mechanical Properties of the damaged beam
Undamaged 10% Damage 25% Damage 40% Damage 50% Damage
Young’s modulus (GPa)
Shear modulus (GPa)
Density (g/cm3)
200 180 150 120 100
143 129 108 86 72
7.84 7.84 7.84 7.84 7.84
1) Apply excitation to undamaged structure.
2) Compute dominant POM of the undamaged structure.
3) Apply excitation to damaged structure.
4) Compute dominant POM of the damaged structure.
5) Compute the difference of the POMs. Damage will be located where the slope of the difference has an evident discontinuity.
Fig. 7. Block diagram describing the steps of the POD damage detection method.
clearest results for detecting the presence of damage in any node. The plots shown in the next sections were obtained by applying the external motion to node 11 of the SSB.
5.3. Different damage positions The effectiveness of the method at locating damage can vary with the position of the damage in the structure. Damage can be introduced in different elements. The effect of this is shown in Figs. 10 and 11, which were generated by introducing a 25% reduction in the mechanical properties in different elements of the mesh. Fig. 10 suggests that if damage is in certain zones of the structure it could be difficult to detect. In particular when damage is located at node 35 of the CB (or closer to the free end) the diagram does not seem to indicate the presence of damage in a clear way. The method seems to work well irrespective of where the damage is located in the SSB.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1356
0.0014
Delta Eigen Vector
0.0009
0.0004
-0.0001 0
5
10
15
20
25
30
35
40
-0.0006
Damage 10% Damage 25% Damage 40% Damage 50%
-0.0011
-0.0016
Element
Fig. 8. Effect of damage severity, CB, damage in element 20.
0.003
Damage 10% 0.0025
Damage 25% Damage 50%
0.002
Delta Eigen Vector
0.0015 0.001 0.0005 0 0
5
10
15
20
25
30
35
40
-0.0005 -0.001 -0.0015 -0.002
Element
Fig. 9. Effect of damage severity, SSB, damage in element 20. Motion imposed to nodes 1 and 41.
5.4. Noise White Gaussian noise was added to previous results of the CB and of the SSB. Figs. 12 and 13 illustrate that as the severity of damage is increased, the tolerance of the newly proposed method to noise also increases. In
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1357
0.00038
-0.00002
0
5
10
15
20
25
30
35
40
-0.00022
-0.00042
-0.00062
Element Damage Node 5 Damage Node 25
Damage Node 10 Damage Node 30
Damage Node 15 Damage Node 35
Damage Node 20
Fig. 10. Twenty-five per cent damage at various locations of the CB.
Damage Node 11
Damage Node 20
Damage Node 31
0.00125
0.00075
Delta Eigen Vector
Delta Eigen Vector
0.00018
0.00025
0
5
10
15
20
25
30
35
-0.00025
-0.00075
Element
Fig. 11. Twenty-five per cent damage at various locations of the SSB, excitation at node 11.
40
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1358 0.0003
0.0002 Effect of noise (10% damage at element 20)
Delta Eigen Vector
0.0001
Noise level 1%
0 0
5
10
15
20
25
30
35
40 Noise level 5%
-0.0001
Noise level 10% Noise level 20%
-0.0002
-0.0003
-0.0004
Element
0.0008 Effect of noise (25% damage at element 20)
0.0006
Noise level 1%
Delta Eigen Vector
0.0004
Noise level 5% 0.0002 0
Noise level 10% 0
5
10
15
20
25
30
35
40
Noise level 20% Noise level 30%
-0.0002 -0.0004 -0.0006 -0.0008
Element
0.0015 Effect of noise (50% damage at element 20)
0.001
Noise level 1%
Delta Eigen Vector
0.0005
Noise level 5% 0
0
5
10
15
20
25
30
35
40
Noise level 10% Noise level 20%
-0.0005 Noise level 30% Noise level 40%
-0.001
Noise level 50%
-0.0015
-0.002 Element
Fig. 12. Effect of noise on damage detection in the CB.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1359
Table 3 Relationship between % value and snr value of the noise level % value
1%
5%
10%
20%
30%
40%
50%
snr
48.417
34.250
28.148
22.046
18.477
15.945
13.980
Figs. 12 and 13 the noise level is indicated by a percent value which is linked to the snr parameter of Matlab according to Table 3. For 10% damage the noise dominates if its percent value is above 10%, in both the CB and the SSB. For 25% damage the tolerance is increased and damage can be seen clearly even at 20% noise level in the SSB, whereas its detection is more problematic in the CB. However, the signal becomes noise dominated for higher levels of noise. Finally, for 50% damage, noise has limited effect on the process of damage detection in both models. While the noise will distort the Mode shapes its effect is small compared to the magnitude of the POM differences.
5.5. Remarks Possibly one of the main drawbacks of the POD damage detection method, as presented in this paper, is that it is based on data coming from a dense grid of sensors distributed on the structure. That would make it expensive and impractical, especially for large structures. The authors believe that the most likely applications of the method under investigation would probably be in the quality control process of industrial production and would allow for a fast and safe detection of structural components (blades, stiffened panels, tubes, etc.) which do not meet the required standards. In such cases it would probably be possible to compute the POMs using data obtained with non-contact optical devices. Another possible field of application could be represented by smart structures already equipped with embedded sensors and actuators which would be used to implement the new technique. Moreover the method could be used also for real-time damage detection, for example in dynamic fatigue tests with no interruption of the test itself [26].
6. Conclusions A new idea for detecting and locating damage in structures has been presented. It is based on concepts of the theory of POD and has been investigated at a numerical level for beams.
Damage is revealed by a sudden change in the slope of the curve showing the value of DPOM, the difference between the values of the dominant POM of the undamaged and damaged structures. An increase in the intensity of damage generates an increase in the value of the change in slope of DPOM. It can be problematic to locate the damage if it affects certain zones of the structure, but in the investigated examples that often happens in the region where stresses and strains are the smallest, and therefore where there is a smaller probability of having damage. The presence of noise can reduce the capability of the method to locate the damage.
The methodology seems promising since it is capable of locating the damage and presents a degree of robustness with respect to the presence of noise. The method is attractive because it is global, easy to implement and can be extended to structures of more complicated shape.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1360 0.0004 0.0003
Effect of noise (10% damage at element 20)
Delta Eigen Vector
0.0002 0.0001
Noise level 1% 0
0
5
10
15
20
25
30
35
40
Noise level 5% Noise level 10%
-0.0001
Noise level 20% -0.0002 -0.0003 -0.0004 Element 0.001
Effect of noise (25% damage at element 20)
Delta Eigen Vector
0.0005
Noise level 1% Noise level 5%
0 0
5
10
15
20
25
30
35
40
Noise level 10% Noise level 20%
-0.0005
-0.001 Element Effect of noise (50% damage at element 20)
0.0024 0.0019
Noise level 1%
Delta Eigen Vector
0.0014
Noise level 5% 0.0009 Noise level 10% 0.0004 -0.0001 0
Noise level 20% 5
10
15
20
25
30
35
40
Noise level 30% Noise level 40%
-0.0006
Noise level 50% -0.0011 -0.0016 -0.0021 Element
Fig. 13. Effect of noise on damage detection in the SSB.
ARTICLE IN PRESS U. Galvanetto, G. Violaris / Mechanical Systems and Signal Processing 21 (2007) 1346–1361
1361
References [1] S. Ness, C.N. Sherlock (Eds.), Nondestructive testing overview, ASTN Nondestructive testing Handbook, Vol. 10, American Society of NDT, 1996. [2] S.W. Doebling, C.R. Farrar, M.B. Prime, A summary review of vibration-based damage identification methods, The Shock and Vibration Digest 30 (2) (1998) 91–105. [3] P. Cawley, R.D. Adams, The location of defects in structures from measurements of natural frequencies, Journal of Strain Analysis 14 (1979) 49–57. [4] O.S. Salawu, Detection of structural damage through changes in frequency: a review, Engineering Structures 19 (9) (1997) 718–723. [5] C.P. Ratcliffe, Damage detection using a modified Laplacian operator on mode shape data, Journal of Sound and Vibration 204 (3) (1997) 505–517. [6] K. Pandey, M. Biswas, M.M. Samman, Damage detection from changes in curvature mode shapes, Journal of Sound and Vibration 145 (1991) 321–332. [7] Y.K. Ho Structural damage identification using changes in vibration characteristics, Ph.D. thesis, Imperial College London, 2000. [8] P. Holmes, J.L. Lumley, G. Bekooz, Turbulence, coherent structures, dynamical systems and symmetry, Cambridge University Press, Cambridge, 1996. [9] J.L. Lumley, Stochastic Tools in Turbulence, Academic Press, New York, 1970. [10] L. Sirovich, Turbulence and dynamics of coherent structures. Part 1: Coherent structures, Quarterly of Applied Mathematics 45 (3) (1987) 561–571. [11] J.P. Cusumano, M.T. Sharkady, B.W. Kimble, Experimental measurements of dimensionality and spatial coherence in the dynamics of a flexible beam, Philosophical Transactions of the Royal Society of London Series A 347 (1994) 421–438. [12] R.V. Kappagantu, B.F. Feeny, Part 2: Proper orthogonal modal modeling of a frictionally excited beam, Nonlinear Dynamics 23 (2000) 1–11. [13] B.F. Feeny, On proper orthogonal co-ordinates as indicators of modal activity, Journal of Sound and Vibration 255 (2002) 805–817. [14] M.F. Azeez, A.F. Vakakis, Proper orthogonal decomposition (POD) of a class of vibroimpact oscillations, Journal of Sound and Vibration 240 (5) (2001) 859–889. [15] X. Ma, A.F. Vakakis, L.A. Bergman, Karhunen-Loe`ve modes of a truss: Transient response reconstruction and experimental verification, AIAA Journal 39 (4) (2001) 687–696. [16] I.T. Georgiou, I. Schwartz, ‘POD (K-L) analysis of the interaction free dynamics of a beam/pendulum configuration: an experimental study’, Proceedings of DETC’01, Pittsburgh, PA, September 9–12, 2001. [17] G. Kerschen, J.C. Golinval, Physical interpretation of the proper orthogonal modes using the singular value decomposition, Journal of Sound and Vibration 249 (2002) 849–865. [18] R. Ruotolo, C. Surace, Using SVD to detect damage in structures with different operational conditions, Journal of Sound and Vibration 226 (1999) 425–439. [19] H. Soon, J.A. Czarnecki, C.R. Farrar, Structural health monitoring using statistical process control, Journal of Structural Engineering (2000) 1356–1363. [20] E. Reissner, On one-dimensional finite strain beam theory: The plane problem, Journal of Applied Mathematics and Physics (ZAMP) 23 (1972) 795–804. [21] P.B. Bornemann, U. Galvanetto, Discrete dynamics of implicit time integration schemes for a dissipative system, Communications in Numerical Methods in Engineering 19 (2003) 737–748. [22] R.R. Craig Jr., Structural Dynamics, Wiley, New York, 1981. [23] P.B. Bornemann, Time integration algorithms for the steady states of dissipative non-linear dynamic systems, Ph.D. thesis, Imperial College London, 2003. [24] G. Violaris, Damage detection in beam structures using the proper orthogonal decomposition, Master Thesis, Department of Aeronautics, The composites Centre, Imperial College London, 2004. [25] D. Hunter, Damage detection in isotropic, cantilever plates using Proper Orthogonal Decomposition, Final Year Project Report, Department of Aeronautics, Imperial College London, 2005. [26] M. Wiercigroch, C.H. Foong, W.F. Deans, An experimental rig to investigate fatigue crack growth under dynamic loading, Meccanica 38 (2003) 19–31.