Crack detection in beam-like structures by nonlinear harmonic identification

The dynamic behavior of beam-like structures with fatigue cracks forced by harmonic excitation is characterized by the appearance of sub and super-harmonics in the response even in presence of cracks with small depth. Since the amplitude of these harmonics depends on the position and the depth of the crack, an identification technique based on such a dependency can be pursued: the main advantage of this method relies on the use of different modes of the structure, each sensitive to the damage position in its peculiar way. In this study the identification method is detailed through numerical examples tested on structures of increasing complexity to evaluate the applicability of the method to engineering applications. The amount of data to obtain a unique solution and the optimal choice of the observed quantities are discussed. Finally, a robustness analysis is carried out for each test case to assess the influence of measuring noise on the damage identification.


INTRODUCTION
he occurrence of cracks in civil or mechanical systems leads to dangerous effects for the structural integrity and causes anomalous behaviors.To not compromise the safety, the detection of a crack in the early stage is of great interest.The presence of a crack not only causes a local variation in the mechanical characteristics of the structure at its location, but it also has a global effect involving the entire structure.For this reason, the dynamic characterization of cracked structures can be used for damage detection in non-destructive tests and, among the various techniques, vibration-based methods provide an effective means of detecting fatigue cracks in structures [1][2][3][4][5][6][7].The vibration-based detection methods exploit the fact that damage in a structure alters its dynamic properties: in particular the presence of cracks can be disclosed by modifications in the linear response (i.e., modal frequencies, damping and shapes associated with each natural frequency, curvature of mode shapes etc.) as well as by the occurrence of nonlinear effects (sub-and super-harmonics, bifurcations, internal resonances etc.).There are two main categories of crack models used in vibration-based detection methods: open crack and breathing crack models.In the first case it is assumed that the crack in a structural member remains open during vibration.This assumption is usually satisfied in notched beams and when the damage is rather large; this model avoids the complexity resulting from nonlinear behavior when a breathing crack is presented.On the contrary, breathing behavior is generally T reported in the case of fatigue cracks, also when the damage affects only a small portion of the cross section of the structural element; it requires a nonlinear model to take into account its effect on the system dynamics.In fact the breathing crack model considers that, during the vibration cycle of a structure, the edges of the crack come into and out of contact, leading to sudden changes in the dynamic response of the structure.Depending on the crack model, vibration based methods are also classified into two categories: the linear and the nonlinear approaches.The first group of methods can identify only the open cracks once at a developed stage, once changes in modal parameters become significant [6,7].For this reason, refined studies focus on the nonlinear response characteristics that can be investigated to identify the presence of the crack in an early stage.In fact, the structures with breathing cracks behave similarly to bilinear systems and hence exhibit nonlinear phenomena in the dynamic response even for low damage.Therefore in the second group of methods the identification can be obtained by assuming as damage indicators some peculiar characteristics of the nonlinear dynamic response such as the presence of sub and super harmonics [8][9][10][11], the changes in the phase diagrams [11,12], the rise of superabundant nonlinear normal modes [13][14][15] and non-smooth bifurcations [16].While these nonlinear effects make the response of beams more difficult to model with respect to notched beams, their appearance clearly marks the boundary between undamaged and damaged behaviour.In fact, as known [17], when a system with a breathing crack is excited by a single harmonic force, distinctive nonlinear features appear in the response.The excitation, in fact, forces the crack to open and close and the resulting clapping of the crack's edges produces harmonics that are integer multiple or fractional multiple of the forcing frequency.These harmonics are commonly referred to as superharmonics and sub-harmonics, respectively.These features are easily detectable when the excitation frequency is in an integer ratio or is a multiple of a resonance frequency of the system; moreover, these would be much more sensitive to cracks characteristics than the modal properties of a linear system.For this reason, increasingly over recent years [11,12,18,19], attentions have been focused on the investigations of the nonlinear effects caused by the presence of breathing cracks and on the associated applications to the problem of damage detection.According to the previous observations, the dynamic behavior of beam-like structures with fatigue cracks forced by harmonic excitation is characterized by the appearance of sub and super-harmonics in the response even in presence of cracks with small depth.Since the amplitude of these harmonics depends on the position and the depth of the crack, an identification technique based on such a dependency has been developed by the authors [19].The main advantage of this method relies on the combined use of different modes of the structure, each sensitive to the damage position depending on the curvature of the mode at its location.In this paper the identification method proposed in [19] is extended and detailed through numerical examples tested on structures with a single breathing crack and of increasing complexity to evaluate the applicability of the method to engineering applications.Dynamic vibrations of small amplitude will be considered so that the crack can be assumed to be stable.The amount of data to obtain a unique solution and the optimal choice of the observed quantities are discussed.Finally, a robustness analysis is carried out for each test case to assess the influence of measuring noise on the damage identification; the robustness of the identification, evaluated through a Monte Carlo simulation, is shown to be quite strong to both measuring and modeling errors envisioning the possibility for in-field applications of this method even in the case of very small cracks.As future developments of this study, another field of investigation will cover the robustness of the procedure with respect to other realistic disturbances such as those caused by constraints imperfections or random distributed mass.Two test cases are studied in the present paper: a simply supported beam and a spanned beam; on the first case, the difficulties in detecting the damage in symmetric systems are addressed, while on the second case the difficulties concerning more complex systems with localized modes are considered.

Generalities
ig. 1a,b present the systems under consideration: a simply supported beam, Fig. 1a, and a spanned beam, Fig. 1b.Both cases have been modeled using standard one dimensional finite elements, while the crack is modeled as explained in the following subsection.It is assumed that only one element has a single sided edge breathing crack and that the switching between open and closed behavior of the crack occurs in the equilibrium configuration.Furthermore it is supposed that the opening and closing of the crack is primarily driven by the flexural dynamics of the beam and that the crack is the only source of nonlinearity in the structure.The considered beams are made of steel with a total length l = 0.6 m, and a constant rectangular cross-section A, of height h=0.04 m and width w=0.02 m; the properties of the material are: E=206 GPa,  =7800 kg/m 3 ,  =0.To characterize the dynamic properties of the structure, distinction should be made between the first natural frequencies, T : . This relation strictly holds for a single-degree-of-freedom oscillator with a bilinear stiffness and is therefore called first bilinear frequency of the system [3].As demonstrated in the literature [3,20,21] the same equation approximates with good accuracy the natural frequencies of a structures with a breathing crack.In conclusion the resonance of the damaged systems occurs at the bilinear frequency i  given by: where  is the i th natural frequency of the undamaged system and d i  is the i th natural frequency of the linear damaged system when crack is always open.

Model of the breathing crack
In several applications, the equilibrium configuration of the system coincides with the switching condition between closed and open behaviour for the crack.Under this condition, the cracked zone of the beam can be modelled with a finite element that has a bilinear element matrix with the discontinuity passing through the origin: this leads to an overall bilinear model where the nonlinear characteristics are independent of the vibration amplitude [3,4,21].The standard element stiffness matrix K for a plane beam element with three degrees of freedom per node, collected in the vector U, once shear deformability is neglected, is described by: where E is the Young's modulus, u I is the moment of inertia of the undamaged beam section, and l is the length of the element.The presence of the breathing crack in the beam is modelled through a particular finite element at the location of the crack.Considering a crack propagating from the upper side of the beam, the damaged element is described by a damaged moment of inertia d I obtained by superimposing to the healthy element another element that has the following stiffness matrix: where  represents the non dimensional flexural damage and ranges between 0 for the healthy section and 1 for the completely damaged section; in Eq. ( 3) the effect of the crack on the axial stiffness is neglected.
Assuming that the opening mechanism is triggered only by the rotations at the nodes of the element as already proposed in [5,18], the breathing mechanism is obtained by the damaged element matrix d K with bilinear behavior: where H is the Heaviside step function dependent on the relative rotations between the sections i and j.In order to compare the results of the proposed model with the results found in literature that often consider as damage parameter the non dimensional crack severity s = a/h, the results are presented by appropriately scaling the damaged axis.The relation between  and s is obtained through a nonlinear compliance function ( ) The particular form of the compliance function   g s is either derived in the literature from theoretical models or from experimental or numerical data.In this work, the compliance function is obtained from a FE model of the open-cracked beam element; it is important to note that ( ) g s depends on the length of the damaged element but the results of the identification will not be affected by the particular choice of the compliance function: this function is merely used to graph the results as a function of s rather than  .
The proposed crack model, while being sufficiently simple, can capture several practical situations and, in particular, fatigue crack are reported to behave in such this way.Moreover, this crack model allows for a consistent reduction of the numerical efforts which is very useful in the solution of the inverse problem.

HARMONIC DAMAGE SURFACES (HDS)
efore tackling the damage identification problem, the direct analysis of the effect of a breathing crack on the dynamics of the beam is developed, by applying the procedure described in [19] and synthetically recalled in the following.First, the amplitude of the sub and super-harmonics of the forcing frequency Ω are evaluated.In practice, as far as a bilinear system is considered, a large harmonic content emerges at a bilinear frequency of the system, Eq. ( 1), when the driving frequency approaches an integer ratio or a multiple of that frequency.The amplitude of the harmonics of interest is obtained by numerically solving the equations of motion of the finite element model of the structure subjected to a point force applied at the loading points shown in Fig. 1.

B
The amplitude of the j-th order super-harmonic is then normalized with respect to the main harmonic amplitude, obtaining the normalized spectrum * ( ) S  around i  .Its maximum value is labeled as ij R and it depends on crack position and severity: * * ( ) max( ), ( / ) The value of ij R depends on the particular choice of the excited harmonic: for instance 22 R is obtained by considering for the second mode the second order harmonic, when the driving frequency is half of the 2 nd frequency.The nondimensional ratio ij R is obviously a function of both the crack depth and its location but, because of the normalization and the bilinear behaviour, is not dependent on the used excitation amplitude, as explained at the beginning of the previous subsection.To characterize a j-th order sub-harmonic it is sufficient to replace in the above equations 1/ j to j ; as an example 2½ R is obtained by considering for the second mode the second order sub-harmonic, when the driving frequency is twice the 2 nd frequency.As the crack position and its severity change, ij R describes a surface that can be computed by performing a parametric analysis.We will refer to this kind of surface as Harmonic Damage Surface (HDS).In particular, Fig. 2a and   /2; a)

Description
or the damage identification purpose, using the information obtained from direct simulations, the excitation and measurement points can be chosen; the test is carried out exciting the system with a sine sweep and measuring only the response signals.The time histories recorded from the sensors are then post processed adopting the same algorithm used for numerical data produced by the FE model to evaluate the corresponding value of m ij R , where the superscript m stands for measured.According to Eq. ( 6), the Fourier transform of the signal is computed and the harmonic peak value around i  is divided by the peak value of the excitation frequency ( / i j  ).
The measured ratio m ij R defines a sectioning plane for the corresponding HDS that intersects it along a curve representing an iso-harmonic-content-curve.This curve collects all the combination of damage-severity and damage-position identified by the performed measure m ij R .As an example, Fig. 2a , where   , ij e p s depends on the particular choice of the HDS and it is a function of both the damage severity and the damage position.In Eq. ( 7), 1, , l n   is the index of the computed n points   , l l p s during the direct analysis to build the corresponding HDS, while 1, , k q   is the index of the q measurements taken on the structure for the identification.Finally, the functional to be minimized represents a global normalized error between the model results and the measured data:

Robustness of the method
In the real world, the identification is affected by measuring errors which originate from different sources such as noise in the measures, environmental excitation, error in the positioning of the sensors etc.The proposed method is intrinsically F robust with respect to random errors, because the m ij R indicators are normalized quantities.In particular, even a large white noise on the measured response will not lead to significant error in the identification.On the contrary, m ij R is sensitive to noise affecting the limited band around either the excitation frequency or the considered super harmonic or other source of nonlinearity present in the system.The amplitude of these errors will be independent of the level of the harmonic ratio due to the damage but they will be dependent on the frequency, on the considered mode, and on the considered harmonic.In [19] the robustness of the method has been tested by considering a cantilever beam and assuming that each measured harmonic m ij R has, independently from the other, an error that is Gaussian, with a mean value equal to zero, and a standard deviation   that is 15% of the average value of the corresponding HDS surface.This level of harmonic error should realistically account for any nonlinearity due to material or non-perfect constraints.The quantification of the error on the identified damage was obtained through a Monte Carlo simulation performed by building a large 50.000samples set of triplets m ij R  , and for each of them the identification was carried out.The robustness of the methodology was evaluated through the statistics on the resulting set of position and severity identified.Since the robustness is also very dependent on the position and on the level of the damage, the Monte Carlo simulation was performed for different combinations of the two parameters.As shown in Fig. 4, it has been found that in all the considered cases the mean value of error is always below 2.5% and reaches 0.1% as soon as the damage severity exceeds 10%; on the contrary, the standard deviation of the error, which measures how the data spread around their mean values, exhibit a different behavior.In the worst scenario i.e. when the damage is far from the constraint and with very low depth (s < 10%), the standard deviation of the position error, reaches 15%, but falls below 1% as the damage severity reaches 10%.The deviation of the identified severity is generally very low, between 1% and 0.5% of the beam height, and is largely independent on the location and severity.The results show that the procedure is able to fully identify damages in any position when s>0.1, and, for very small damage (0.05<s<0.1), the position remains uncertain but the severity is still well assessed.Moreover, results show that in general, the dispersion of all the 50000 data samples of the Monte Carlo simulation remains quite close to the correct value improving confidence on the proposed technique.

Generalities
he present study deals with specific aspects and characteristics of the identification method.In particular, three different numerical tests are carried out to highlight each a specific problem: a) improving the identification by using more than one sensor; b) how to identify the damage in a symmetric structure; c) how the identification T changes as the structure becomes more complex.It is here important to recall that the method detects the damage level through the response of a selected mode of the systems thus, by considering just one mode, a large damage close to a curvature node of that mode has the same effect of a small damage in a curvature antinode.

Damage detection in symmetric structures
When a symmetric structure is considered, in this case a simply supported beam, the are either symmetric or anti symmetric.Thus, in principle, it is difficult for the harmonic identification method to distinguish between the real damage position and its symmetric, while the correct evaluation of the damage severity should not be affected.In practice, however, the presence of the damage breaks the symmetry of the system and this can be exploited in the identification by using a non symmetric excitation of the structure.In particular in Fig. 5

Information from sensors located at different points
The effect of the position of the sensors is illustrated using the simply supported beam.Fig. 7a-g shows the results of the identification for different positions of the sensors p s =[0.12 0.26 0.32 0.4 0.6 0.68 0.8] the vertical line shows for each plot the sensor position.As it can be seen all the sensors highlight the correct position and severity of the damage, however some of them, which in general depend on the relative position between damage, force, considered mode and position of the sensor, can produce some other zone of damage.For example, in the particular case of Fig. 7e, with the sensor located at p s =0.6, a spurious zone of detected damage around p=0.85, s=0.2 appears.In order to avoid this misleading result, the use of an average between the several results provided by all the sensors can be used, Fig. 7h.respectively.It is clear that more information is available at a reasonable low cost being, in general, quite inexpensive to use more than one sensor in the analysis.Moreover, the use of 8 sensors allows an identification that is able to locate the damage, unless it is close to a constraint, from 10% severity down to just 5%.In this paper, the problem of symmetric structures is also addressed and it is shown how, even the small asymmetry caused by the presence of the damage can be exploited, by exciting the system with a non-symmetric load to correctly obtain the damage identification.Moreover, the use of more than one sensor on the structure is investigated, considering in the minimization process the information from each measured point.The obtained results, show that with sensors in different positions the method is able to locate the damage, unless it is close to a constraint, starting at just 5% severity Finally, the case of a more complex structure with multiple constraints is addressed.The tests show that using three modes is not sufficient to identify the damage at any location.The use of a fourth mode is necessary to identify the damage in any position.This aspect will be the subject of further investigation, increasing also the complexity of the system; moreover, it is planned to verify with an experimental campaign the effectiveness of the proposed method on the controlled laboratory test cases.
3. The beams have a breathing crack of depth a at distance d from the left-end constraint; p = d/l is the non dimensional position

1  1 u T and 1 d
two constituent sub-models (open crack, superscript d, and closed crack, superscript u)and the first natural frequency of the system: this frequency can be obtained as the average of the natural periods of the two sub-models

Figure 2 :Figure 3 :
Figure 3: Spanned beam: Harmonic damage surfaces R32 order 2 super-harmonic component for Ω = 3 reports a representation (dashed line) for the HDS surface of 22 R furnished by the model sectioned by the plane 22 m R = 0.0135; this value corresponds to a damage located at p = 0.27 and with severity s = 0.1, that is a point of the curve..The damage identification is obtained by computing, the error function ij e for each ij R surface p and s for which ( , ) f p s is minimum, represent the position and the severity of the damage.

Figure 4 :
Figure 4: Standard deviation of the errors on the identified damage as the damage position and the severity changes for a 15% Gaussian random error on the measured harmonics.a) Standard deviation of the position error; b) Standard deviation of the severity error.

Figure 7 :
Figure 7: Nonlinear harmonic identification on a simply supported beam with damage p = 0.27, s = 0.1 excited at l F /l=0.06.a)-g) Results with one sensor located at different points; h) average of the results provided by all the cases a)-g).In order to account for several input sensors the

Figure 8 :
Figure 8: Nonlinear harmonic identification on a spanned beam with damage p=0.4,s=0.1 excited at l F /l=0.5.a) Three mode identification; b) Four mode identification.