Microsoft Word - numero_29_art_27 P. Casini et alii, Frattura ed Integrità Strutturale, 29 (2014) 313-324; DOI: 10.3221/IGF-ESIS.29.27 313 Focussed on: Computational Mechanics and Mechanics of Materials in Italy Crack detection in beam-like structures by nonlinear harmonic identification Paolo Casini, Fabrizio Vestroni University of Rome la Sapienza p.casini@uniroma1.it, vestroni@uniroma1.it Oliviero Giannini University Niccolò Cusano oliviero.giannini@unicusano.it ABSTRACT. 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. KEYWORDS. Fatigue crack; Damage identification; Nonlinear dynamics; Cracked beam model. 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-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 P. Casini et alii, Frattura ed Integrità Strutturale, 29 (2014) 313-324; DOI: 10.3221/IGF-ESIS.29.27 314 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-11], the changes in the phase diagrams [11, 12], the rise of superabundant nonlinear normal modes [13-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 super- harmonics 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. SYSTEMS UNDER INVESTIGATION 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/m3,  =0.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 F P. Casini et alii, Frattura ed Integrità Strutturale, 29 (2014) 313-324; DOI: 10.3221/IGF-ESIS.29.27 315 and s = a/h is the non dimensional severity of the crack. The simple supported beam, Fig. 1a, has been forced under two different loading condition: almost symmetric 0.45F l l      and asymmetric 0.06F l l      . Also for the spanned beam two different locations for the load are considered: 2 3 Fl l and 2 5 Fl l , where 1 2l l l  . The beams are modeled using standard one dimensional Euler-beam elements, with three degrees of freedom for each node, as well as the damaged element; the mesh consists of 60 finite elements. a) b) Figure 1: Systems under investigation. a) Simple supported beam, l=0.6 m, lF=0.26 m (quasi symmetric) or lF=0.04 m (asymmetric). b) Spanned beam, l1=0.4 m, l2=0.2 m, lF=0.3 m or lF=0.5 m. To characterize the dynamic properties of the structure, distinction should be made between the first natural frequencies, 1 1 2u uT    and 1 1 2d dT    , of the two constituent sub-models (open crack, superscript d, and closed crack, superscript u) and the first natural frequency 1 of the system: this frequency can be obtained as the average of the natural periods of the two sub-models 1 uT and 1 dT : 1 11 1 1 1 4 u d u d u d iT 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: u d i i i u d i i        (1) where ui is the ith natural frequency of the undamaged system and d i is the ith 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: 2 2 2 / 0 0 0 0 /  6 / 0  12 / 12 / 6 /     4     0 6 /     0 ,             / 0     0       12 / 6 /              4 iu u i i u ju j j uA I A I vll l l lEI K U uA Il vsymm l l                                   (2) P. Casini et alii, Frattura ed Integrità Strutturale, 29 (2014) 313-324; DOI: 10.3221/IGF-ESIS.29.27 316 where E is the Young's modulus, uI 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 dI obtained by superimposing to the healthy element another element that has the following stiffness matrix: 2 2 2 0 0 0 0 0 0     6 / 0  12 /    12 / 6 /           4   0 6 /     0 ,               0 0     0                   12 / 6 /              4 u u d c u ll l l lEI I I K l I symm l l                       (3) 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 dK with bilinear behavior:     1,              0,   i j d i j c i j i j K K H K H                 (4) 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 ( )g s : ( )u d u I I g s I    (5) 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 P. Casini et alii, Frattura ed Integrità Strutturale, 29 (2014) 313-324; DOI: 10.3221/IGF-ESIS.29.27 317 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 ijR and it depends on crack position and severity: * * ( ) max( ),       ( / ) i ij i S R S S S j     (6) The value of ijR depends on the particular choice of the excited harmonic: for instance 22R is obtained by considering for the second mode the second order harmonic, when the driving frequency is half of the 2nd frequency. The non- dimensional ratio ijR 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 2nd frequency. As the crack position and its severity change, ijR 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 2b present the HDS for 22R for the simple supported beam obtained for symmetric and asymmetric loading condition respectively; it is worth noticing that in the first cases the ijR coefficients are considerably smaller, since the loading position is close to the modal shape node. a) b) Figure 2: Simple supported beam: Harmonic damage surfaces R22 order 2 super-harmonic component for = /2; a) symmetric loading condition; b) asymmetric loading condition. Similarly, Fig. 3a and 3b present the HDS for 32R for the spanned beam obtained for two different load locations: 0.3Fl  m and 0.5Fl  m. It can be observed that the loading conditions can significantly affect the HDS. The amplitude of ijR reflects the influence that the damage has on the particular mode: therefore, the surface has a monotonic behaviour with respect to the damage severity whereas, along the p axis, the surface should resemble the curvature of the corresponding mode. In this paper the HDS related to 12R , 22R , 32R are used for the simple supported beam and 12R , 22R , 32R , 42R for the spanned one. P. Casini et alii, Frattura ed Integrità Strutturale, 29 (2014) 313-324; DOI: 10.3221/IGF-ESIS.29.27 318 a) b) Figure 3: Spanned beam: Harmonic damage surfaces R32 order 2 super-harmonic component for Ω = 3 /2; a) 2 3 Fl l ; b) 2 . 5 Fl l NONLINEAR HARMONIC IDENTIFICATION METHOD 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 mijR , 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 mijR 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 mijR . As an example, Fig. 2a reports a representation (dashed line) for the HDS surface of 22R furnished by the model sectioned by the plane 22 mR = 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 ije for each ijR surface       ,1 1,,, ,      n q mm ij l l ij kij ij m l k ij ijm ij R p s RR p s R e p s R n qR         (7) where  ,ije 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 lp 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: 2 , ( , ) ( , )ij i j f p s e p s  (8) the values of p and s for which ( , )f p s is minimum, represent the position and the severity of the damage. 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 P. Casini et alii, Frattura ed Integrità Strutturale, 29 (2014) 313-324; DOI: 10.3221/IGF-ESIS.29.27 319 robust with respect to random errors, because the mijR 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, mijR 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 mijR 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.000 samples set of triplets mijR , 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. a) b) 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. 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