Modal curvature-based damage localization in weakly damaged continuous beams Mechanical Systems and Signal Processing

Modal curvatures have been claimed to contain local information on damage and to be less sensitive to environmental variables than natural frequencies. However, simply using the difference between modal curvatures in the undamaged and damaged states can result into localization errors, due to the complex pattern that this quantity presents when considering broad damages or higher order modes. In this paper, we consider weakly damaged con- tinuous beam and we exploit a perturbative solution of the beam equation of motion to obtain an analytical expression of the modal curvature variations in terms of damage dis- tribution. The solution is then used to introduce a ﬁltering technique of the modal curvature variation and to set up the inverse problem of damage localization based on modal curvatures only. Using numerical examples and experimental tests, we show that modal curvatures can be used for precise damage localization, once properly ﬁltered. (cid:1) 2018 Elsevier Ltd. All rights reserved.


Introduction
Structural health monitoring techniques based on the measurement of modal response have been the subject of much research in the last decades [1][2][3][4] fostered by the availability of low-cost, reliable sensors suitable to monitor large civil constructions including buildings, bridges and aqueducts [5][6][7][8][9][10].
The change of natural frequencies has been one of the first approaches used in damage detection [11][12][13][14][15][16], thanks to the easiness and robustness of their measurement in comparison to other modal quantities. The intrinsic drawback in the use of modal frequencies is, though, their limited sensitivity to local damage of moderate severity, which might lead to significant errors in the identified damage parameters. This can be particularly critical in the presence of environmental factors like variations of temperature or humidity, which can cause changes of natural frequencies of the same order of magnitude of those induced by damage [17,10]; as such, frequency changes can be erroneously related to the presence of damage. What is more, the inverse problem of damage detection techniques based on eigenfrequencies is often ill-conditioned [2].
Some of these limitations can be overcome by using the changes of modal displacements [18] and curvatures [19,20], which are claimed to be more sensitive to local damage. Such an approach was explored in the last decade of the last century by Pandey et al. [21] and, since then, has given rise to quite a number of follower studies investigating the use of modal curvatures as to solve the problem of damage localization and assessment in beams [19,[22][23][24][25][26], arches [27], and beam-like structures [28][29][30][31]. Compared to other localization techniques, modal curvatures have several advantages. In fact, for narrow damage, such as a notch, the variation of modal curvatures localizes in the neighborhood of the damage [21,27]. Moreover, curvatures are less sensitive to environmental factors than natural frequencies. These favorable properties lead us to examine more in-depth the effectiveness of modal curvature variations and to extend a procedure for damage localization initially proposed in [20] for single span beams, and based on the change of modal curvatures between the undamaged and damaged states. Although appealing, the use of curvatures for damage localization has to be pursued carefully. To begin with, an important point is their actual measurement. The derivation of modal curvatures from modal displacements by means of numerical differentiation is an ill-conditioned problem. In the literature, several procedures have been presented to deal with this issue, including spline interpolation [32], wavelet transform [33], modified Laplace operator [34] or space-wavenumber Fourier transform [35]. All these techniques require a high number of measurement points to obtain reliable values of modal curvatures and to be accurate in damage localization [29,33]. A direct measurement of curvatures can be substantially obtained through the time-history of strains on the surface of the beam, which may be related to curvatures by accepting the hypothesis that cross-sections remain plane. However, in this case, due to the expression of the frequency response function, normalized mode shapes can be obtained only with some approximation [29,36].
In this paper, a perturbative solution of the dynamics of the damaged Euler-Bernoulli beam is used to determine an analytical approximation of the relationship between the difference of modal curvatures in the undamaged and damaged states and the damage distribution. In this way, the obtained closed form expression highlights the different contributions that make the modal curvature difference not concentrated at the damage location as, indeed, one would expect [21,37], and point out the necessity of the modal curvatures to be filtered by these spurious terms. Several damage cases are considered including single and multiple damages as well as narrow and extended. The applicability of the proposed procedure is demonstrated with reference to experimental tests carried out by the authors.

A perturbative solution for the dynamics of continuous beams
An analytical expression of the transverse modal response of a weakly damaged continuous beam governed by dynamic Euler-Bernoulli equation is derived by a perturbative approach. A straightforward expansion of the solution is carried out in terms of the smallness parameter e that represents the damage intensity.
The governing equations of the transverse motion of the n-spans beam shown in Fig. 1 in terms of the parameters a 4 i ¼ q i ' 2 i =EI i , and of the damage function g s ð Þ (max s2 0;1 ð Þ g s ð Þ ¼ 1), representing a stiffness variation, here assumed localised in the i-th span. It is also assumed that g s ð Þ is a smooth function with compact support, i.e., it is non-zero only in the damaged region. Here and henceforth, an overimposed prime will denote differentiation with respect to the dimensionless abscissa s, whereas an overimposed dot will indicate the time derivative.
The boundary conditions are prescribed to guarantee continuity of displacements/rotations and balance of forces/moments at each internal node whereas at the two ends of the beam hinged support or free boundary conditions are prescribed ( Fig. 1 In view of the positiveness of the coefficient of the leading order terms, the solution of the PDE system (1) can be obtained by separation of variables; namely the transverse displacement of the i-th span can be written as which upon substitution in (1), by simple passages, gives the following system of ODEs which is written for the sake of conciseness in matrix form, by defining: / s ð Þ ¼ / 1 s ð Þ; . . . ; / n s ð Þ f g the mode shape, x the vibration frequency, I the n Â n identity matrix and K and N two diagonal matrices given by One normally looks for a solution of (4) in the form that represents the deflected shape of the i-th span. On substituting (6) into (2), one obtains a homogeneous system of equations in the 4Ân unknowns a 1 ; b 1 ; c 1 ; d 1 ; . . . :; a n ; b n ; c n ; d n f g whose non-trivial solution is found if and only if the determinant of the coefficient matrix vanishes. This condition gives a transcendental equation in x 2 whose infinite solutions represent the modal frequencies of the continuous system. Once these characteristic frequencies are computed, one can rewrite (4) in terms of the k-th mode shape / k and modal frequency x k : In the following developments, we assume that both the mode shape / k and the frequency x k are analytic function of the damage intensity e and that x k is a simple eigenvalue. Following a classical procedure of perturbation techniques, one can carry out a straightforward expansion for small damage intensity (see for instance [38,20]): which, upon substitution into (4), lead to the following cascade of ODEs together with the boundary conditions (2) expanded in series as well. Eqs. (9) and (10) follow the classical structure of a multiscale differential problem: the zero-th order Eq. (9) represents the modal response of the undamaged beam; the first order Eq. (10), is governed by the same differential operator of the intact system, but has a forcing term dependent on the zero-th order solution e / k . The solution of first order Eq. (10) / k can be expressed as sum of a homogeneous solution, that still satisfies the boundary conditions of the zero-th order problem, and a particular solution. On passing we also note that the mode shape vectors of the intact system e / k s ð Þ satisfy the orthogonality condition (see A) and that the differential operator in (9) is self-adjoint with respect to the scalar product Á; Á h i just introduced; thus the eigenvectors are an orthogonal basis in the space of continuous and differentiable vector functions with support in 0; 1 ð Þ. With this in hand, the solution of the first order problem can be expanded in terms of the eigenvector basis e / 1 ; e / 2 ; . . .
where the coefficients of the series b kl define the contribution of the l-th mode of the undamaged structure to the first order variation of the k-th mode. The expression of / k in Eq. (12) leads to a norm of the first order approximation / k equal to 1 À e.
In addition, it is a matter of simple calculations to show that (12) satisfies the boundary conditions (2). In view of (12), the solution of the differential problems (9) and (10) can be computed once the coefficients b kl of the series and the first order eiqenfrequency variation x k are known. The following derivations are devoted to calculate these two quantities.
We start by substituting Eq. (12) into (10) and noting that Eq. (9) holds for e / k , we obtain: By multiplying the previous equation by e / T j s ð Þ and integrating between 0 and 1, we get having taken into account the mode orthogonality condition (11). The integral on the right hand side of the equation is rewritten as where the boundary term vanishes because the damage shape function g s ð Þ and its derivatives vanish at the boundaries. By using such a result together with Eq. (9) written for e / l , we arrive at b kj e which gives, for j ¼ k, the coefficients of the series (12). Eq. (17) expresses the requirement that the zero-th order solution must be orthogonal to the forcing term of the first order equation, i.e., right-hand side of (10), which is a singular equation for / k ; this requirement indeed represents the Fredholm solvability condition for the first order problem (see [38]).
Eqs. (17) and (18) give the first order perturbative solution of the direct problem for an assigned damage function N s ð Þ, which allows us to obtain an analytical expression of the modal curvature variations due to damage.

Damage localization through modal curvatures
The inverse problem of damage localization based on the experimental measurements of the modal curvatures is introduced in this section. It is shown that one can determine the k-th modal damage shape, that is the damage shape multiplied by the k-th modal curvature, and not the damage shape itself, that must be reconstructed from the modal contributions.
In the following, we will indicate the k-th modal curvature, i.e., the second derivative of the k-th mode shape, as w k s ð Þ :¼ / 00 k s ð Þ, whereas the k-th modal damage shape as where we have followed the notation of previous section for the modal curvature perturbation w k s  (2); accordingly the modal curvatures satisfy the orthogonality condition: that is an immediate consequence of the modes orthogonality. Therefore one can compute the norm of the modal curvatures of the undamaged beam as: Eqs. (17) and (18) can then be written in terms of modal curvatures only: From the experimental perspective, it is particularly advantageous to deal with normalised quantities. Thus, by using (21), we can compute that in view of (12) allows us to write the k-th normalised modal curvature difference as Indeed by exploiting the fact that modal curvatures are an orthonormal basis, we can write the function N s ð Þ e w k s ð Þ as that represents a projection of the damage function over the k-th curvature. Eq. (26) can be multiplied by e w T i ; ik, and integrated between 0 and 1 to give in which we made use of the norm (21) and of the definition (18). By using the definition (19) together with (27), we can write (26) as Due to (25) where Dx 2 k ¼ e k is k-th modal frequency variation, between the undamaged and damaged systems. Eq. (29) highlights that the modal damage shape is related to the modal curvature variation, the measured quantity, but several other contributions are present. In general, the modal curvature variation is sensitive to the damage position, the same way as for the frequency variations; in particular, if the damage occurs in a region where the k-th modal curvature vanishes, the k-th modal curvature variation does not give any information on the damage shape. In addition, Eq. (29) shows that the damage shape can be evaluated only after the k-th modal curvature difference is filtered as that equation suggests; interestingly, all the quantities that appear on the right hand side of the equation are indeed easy to measure: Dw k is the difference between the normalised k-th modal curvatures of the damaged and undamaged beam, Dx 2 k is the difference of the corresponding modal frequencies, e w l is the l-th undamaged modal curvature divided by its norm k e w l k. Dealing with normalised quantities has the further advantage that one does not need to calculate the modal mass of the system which is often difficult to assess.
In order to evaluate the different terms that cause the modal curvature difference to have multiple peaks outside the damage region, one can apply the same technicalities to invert Eq. (29) and obtain the expression of Dw k in terms of g k .
By making use of (25), we can write which again highlights the fact that apart from the damage shape g k , the k-th modal curvature difference contains the contributions of all modal curvatures. In particular, by writing Eq. (30) for the q-th span in which no damage occurs, one obtains which once more shows that the modal curvature difference in the q-span is non-zero despite the damage being located at the i-th span. This result further highlights our statement that, without proper filtering, the modal curvature difference is not an effective quantity to be directly used for damage localization. The accuracy in evaluation of the modal damage shape depends on the accuracy in evaluating the three terms of the right hand side of Eq. (29), which is mainly related to the available eigenfunctions to be used in the summation (third term). Before applying the introduced filtering technique to experimental data on a continuous beam, the simple pseudoexperimental problem of a free-free beam is considered. For this problem most of the quantities in Eq. (29) can be obtained in closed form, avoiding the additional numerical errors of computing the modal curvatures from the modal displacements, thus focusing on the number of modes in the summation in Eq. (29).
In Fig. 2 three damage cases are considered whose shapes are represented by the grey area and whose intensity was e ¼ 0:05 for all damages. For each case, selected modal damage shapes are reported for an increasing number of modes in the summation (coloured thin lines), together with the absolute value of the unfiltered modal curvature difference (thick blue line). In the first damage case considered, the direct use of the variation of modal curvature still gives a satisfactory result, whereas, in all the other cases, significant peaks arise in the undamaged area, whose values are comparable with those occurring at the correct location of damage. In most cases, the filtering technique yet reduces these spurious peaks leading to a unique definition of the modal damage shape and a unique localization. In particular, for the case of two damages (Fig. 2  (b)), the first modal curvature variation has a significant peak in the undamaged area between the two damages that the filtering technique eliminates. The filtered second and third modal damage shapes are shown in Fig. 2(c) and (d) in the case of a broader damage. For case (c), the observation of g 2 would be enough to localise the damage, yet g 3 is reported to stress that the modal damage shape is zero where the corresponding modal curvature (the third in this case) vanishes, leading to a double peak in the filtered quantity. This also causes that, whatever number of modes is used to reconstruct g 3 , the correct location of the damage, i.e., the maximum in the modal damage shape, is slightly shifted with respect to the actual damage location. All these circumstances suggest to introduce a damage index which combines the available modal damage shapes, as proposed in the next section.

Experiments
The proposed procedure was experimentally verified by testing ABS 3D-printed beams at the laboratory of the Department of Structural and Geotechnical Engineering, Sapienza University of Rome. Four identical specimen were 3D-printed with different notches to simulate different damaged configurations (see the details in Table 1).
The geometry of the beam, as well as the material properties are displayed in Fig. 3, whereas pictures of the actual specimen and a detail of the damage are shown in Fig. 4(a) and (b). The beam was suspended with a soft rubber band and connected to an electrodynamic shaker for excitation by a metallic stinger with a circular cross-section 1 mm in diameter (Fig. 4  (c)). The ratio between the beam and the stinger bending stiffnesses was about 100 and, hence, the stinger acted like a hinge, meaning that the rotational constraint practically vanishes and the specimen can be modelled as a continuous beam. The shaker was driven by a white noise input signal up to a frequency of 2500/3000 Hz. To measure the proper frequency response function, the input force was recorded with a load cell installed at the base of the stinger. The response was measured using seven strain gauges equally spaced along the beam axis, in the damaged and undamaged cases. The beam was slender enough (h=L < 0:1, with h height of the cross section and L beam span), to relate the strain e i measured by the i-th strain gauge at the top flange to the curvature of the center line by w i ¼ 2 e i =h. The sensors TML FLAB-6-11 foil resistance strain gauges, with a resistance of 120 Ohm and maximum strain 5%, were installed in a quarter-bridge circuit. The sensitive pattern of the strain gauge was 6 mm long, which permits to recover a precise measurement of local strain. This gauge length has a frequency response limit of around 250 kHz; the data were recorded using an acquisition system with a 16-bit A/D converter, at a sampling frequency f s ¼ 25600 Hz.
The experimental frequency response function of the system was obtained by computing the cross-power spectral density between the measured input and output signals with a resolution of 0.78 Hz, corresponding to a windowing length of 1.28 s. The frequency response function in terms of strains is tied to the modal parameters as follows: where f k is the damping of the k-th mode. Assuming that damage does not affect mass distribution, q is constant, and Þ 2 dx ¼ qLjj/ k jj 2 represents the modal mass. Differently from what happens when the response is measured in terms of accelerations, only an approximate mass normalization is obtainable if strains are measured [29,36]. However, the use of modal curvatures normalized to unit unbinds us from these approximations.   To summarize, the proposed damage identification methodology relies on the following steps: 1. the time-history of the modal strains is measured through strain gauges by exciting the structure with a broad spectrum white-noise; 2. the frequency response function is calculated and the modal curvatures and modal frequencies are extracted from it; 3. the modal curvature difference Dw k is computed and used into Eq. (29) to obtain the damage shape function g k .
The different damage cases are described in Table 1, with s location of damage and h D residual height (see Fig. 3). Cases D1 and D2 include damage located just under the strain gauge #2 with two different depths of the notch, respectively with residual height h D ¼ 3h=4 and h D ¼ h=2 of the undamaged height. Case D3 includes one damage located in between strain gauges #5 and #6 with h D ¼ h=2. Case D4 has two notches, one located under strain gauge #2, plus a second notch located between strain gauges #5 and #6, both with h D ¼ h=2. The values of the first four identified natural frequencies are reported in Table 2 for the undamaged case (U), and the four cases of damage. Notice that a reduction of natural frequencies with increasing damage is observed, but for mode 2, with the slightest damage D1. Such an increase is sometimes observed for a low stiffness reduction, which inevitably accompanied by a reduction of mass. The combination of these two effects, along with the experimental errors, may conceal the reduction of frequency. The reduction of the natural frequencies is also shown by the shift of the second peak in Fig. 4 (d), which compares the modulus of the frequency response function in terms of strains in the undamaged and damaged (D2) cases. Fig. 5 shows the normalized modal curvatures for the first four identified modes, in the undamaged and damaged states. Also in this case, the increase in modal curvature is rather evident when the damage is located under the strain gauge, while it is less prominent when damage is located in-between two strain gauges (cases D3 and D4). However, the curvature not only increases locally in the neighbourhood of damage, but presents a complex pattern, especially for higher order modes. This is shown in the left column of Fig. 6, which reports the differences between modal curvatures in the undamaged and damaged states for the four identified modes. In fact, modes 2 and 4 for the case D1 present peaks under strain gauge #5, where there is no damage, moreover, mode 4 for the case D2 has a peak under strain gauge #6, where again no damage is present, together with a shift of the first peak. What is more, when damage is located in between two strain gauges (cases D3 and D4), the change of curvature tends to peak in a very large area between strain gauges #4-6, so that a precise damage localization becomes difficult; this pattern may become even more involved for broad damages as shown in [20,25]. The modal damage shape, computed through Eq. (29), is reported on the right column of Fig. 6 for the first four modes. When compared with the left column, the peaks of the modal damage shape exhibit smaller amplitudes outside the damaged area and, in general, present narrower peaks.
In order to use the information available from all modes and to reconstruct the damage shape g s ð Þ, we introduce an overall index of damage, that we call Damage Index (DI): that at position k weighs the modal damage shape multiplied by the inverse mode number. This approach privileges the information from low order modes that are, in general, more reliable. A similar function was used in [39,40] for the unfiltered modal curvatures. In Fig. 7, the results obtained by applying Eq. (33) to the experimental data are compared with other localization techniques: a Pandey-like approach based on the application of Eq. (33) to the unfiltered modal curvature differences (second row in the array plot), an energy-based approach proposed in [41] which we call Energy Index (third row) and the Coordinate Modal Assurance Criterion (COMAC) (see for instance Ref. [42]) (fourth row). All indices have been rescaled between 0 and 1 to properly compare them; as such, the identified location of damage is at the position where the value of the index is closer to 1.
In accordance to the results in Fig. 6, the proposed approach gives a clear indication of the damage position in all cases. For damage D1, which is the weakest damage examined, the procedure gives an accurate localization of the damage, whereas all other techniques give a false indication of the damage in other position, most probably due to the peaks in the second and fourth unfiltered modal curvature differences seen in Fig. 6. This problem is partially overcome in the case D2 of a more intense damage for all the identification techniques but COMAC being non-zero in an area much wider than the actual damage width. When the damage is in-between two strain gauges (case D3), both the Damage Index and the modal curvature difference have a maximum in correspondence of the strain gauge 5; however the Damage Index has the further advantage of being more concentrated around the actual damage location, whilst the modal curvature difference predicts a wider damaged area. In the case of two damages (D4), the proposed index locates both the damaged region whereas all the other techniques failed.

Conclusions
Curvature variation is a quantity that has recently gained interest in damage detection procedures for its local information content and its limited sensitivity to environmental effects. Nowadays, thanks to the availability of low-cost sensors, modal curvatures can be measured directly. However, each modal curvature has its own sensitivity to damage location, the same way as for modal frequencies. To clarify some of the theoretical problems underlying the use of modal curvatures in damage localization, we have obtained a perturbative solution of the equation of motion of a continuous damaged Euler-Bernoulli beam. With this in hand, we have derived an explicit relationship between modal curvature variation and modal damage shape. The latter is the damage function multiplied by the modal curvature shape and is apparently related to the variation of modal curvatures, yet it contains other terms that could play a significant role. This is particularly important in damage localization procedures based on curvatures only, as significant localization errors can occur, if the modal curvature variation is not properly filtered. The relevance of these spurious terms depends mainly on the extension of the damage, and its location with respect to the position of measurement points. The former is assessed by considering a pseudoexperimental case on a free-free beam, that can be fully treated analytically, thus avoiding the additional numerical errors introduced by the discretization. For what concerns the latter, experiments have been carried out on a beam with different damage configurations. Several localization techniques have been compared and their capabilities in terms of damage location were assessed.
In conclusion, when the damage is narrow and the modal curvatures are measured close to the damage, as done in [21,37], the modal curvature variation gives accurate enough localization of the damage. In all other cases, modal curvature variations have to be properly filtered by applying the proposed technique.