Analysis and Damage Identification of a Moderately Thick Cracked Beam Using an Interdependent Locking-Free Element

The Timoshenko interdependent interpolation element, based on the assumption of cubic interpolation for the transverse displacement and quadratic interpolation for the rotation, is developed for both the static and the dynamic problems. Next, the different behavior of a beam due to the presence of a damaged zone is investigated and the problem of identifying diffused crack affecting a portion of the beam using natural frequencies is studied. The damaged zone can be completely taken into account by introducing only three parameters, and for the inverse problem, numerical optimization is applied to define their values.


Introduction
There now exist a large number of publications devoted to the task of identifying diffused crack affecting a beam [1,2]. The increasing attention about this problem in the last decades is essentially due to the growing number of applications in structural University of Pavia, Pavia, Italy engineering such as the study of the stability of existing structures hit by earthquake or the monitoring of structural integrity [3].
This paper focuses on a diffused crack identification method based on the natural frequencies that uses the Timoshenko interdependent interpolation element [4][5][6][7]. The main peculiarity of this element is its capability to avoid the shear locking problem, which may affect discretization approaches for moderately/thick beams that are hereinafter intended as shear deformable beams in the spirit of Timoshenko beam theory. In particular, the interdependent interpolation element is obtained by using the Hermite cubic functions for the transverse displacement and interdependent quadratic interpolation for the rotation. By doing so, the element exhibits a locking-free behavior also in the thin beam limit. This feature is not shared by the linear Timoshenko beam element that is well known to be affected by the locking phenomenon when the beam depth h tends to zero.
Structural health monitoring is nowadays indispensable for the reliable assessment of structures under both service loads and some extreme loads such as earthquakes and strong winds. Although recently many numerical techniques such as changes in modal data [8][9][10], the transfer function parameters [11], the flexibility matrix [12], the strain energy [13], the frequency response functions [14,15] and perturbation approaches [16,17] have been developed to localize and to define the magnitude of diffused cracks affecting structural elements, many of them continue to require experimental modal data of both the damaged and the undamaged elements.
The presented damage identification method, instead, is based only on the comparison between the natural frequencies resulting from experimental studies and those obtained by investigating the dynamic behavior of the same undamaged beam through numerical computation, i.e., on the definition of the differences between the real structure and the corresponding design model. Under the hypothesis that the initial state is fully known (and undamaged), the axioms of structural health monitoring set in [18] are met and referring specifically to axiom 2 the two systems under comparison are the real structure and the one simulated numerically via the approximation approach presented herein.
As suggested in many papers [19][20][21][22][23][24][25], this is relatively easy to do for various reasons: the simplicity of measuring natural frequencies (mode shapes requiring a very large number of sensors and affecting by measurement errors) and, under the assumption of linear behavior, the opportunity to describe the diffused crack affecting a beam with only three parameters. These three parameters are the position, the extension and the coefficient of stiffness reduction (or magnitude), and very few frequencies are required to define their values. Therefore, given this scenario and well aware of the huge amount of methods available for identification of structural damage, the ultimate goal of this paper is to provide a damage identification method for moderately thick beams, which enjoys the following distinctive features: 1. the method is numerically inexpensive in that it does not call for a massive amount of data to be initialized and work properly; 2. the method is based on a relatively simple damage model that does not call for damage mechanics or nonlinearity whatsoever. It should in fact be kept in mind that the overall goal is to determine damage amplitude and location, whereas no attention is paid to post-damage response that would call for finite-step nonlinear analysis or the like; 3. the quality of data needed for the proposed method to converge (natural frequencies) is well known to be the most easily available in the field of structural vibrations, and this should be considered a dramatic pro. 4. Last but not least, the adoption of quite a sound finite element able to rule out locking pathologies from the formulation should be considered a distinctive feature of the proposed approach.

Interdependent Interpolation Element: Formulation
As well known, the Timoshenko beam element is a beam that presents shear deformation. So plane sections normal to beam axis remain plane, but not necessarily normal to the axis during the deformation. According to the classical stress theory, the strain energy U for an isotropic linear elastic material occupying region Ω can be written as where ε is the strain tensor and σ is the Cauchy stress tensor conjugated to ε. It is expedient to remember that the strain tensor must satisfy the compatibility equation and that the constitutive equation is given by where λ and μ are the Lamé's constants that uniquely determine the material isotropic behavior.
Assuming an orthogonal Cartesian coordinate system (x, y, z), in which the x-axis is coincident with the centroidal axis of the undeformed beam, the y-axis is the neutral axis and the z-axis is the third axis, the displacements field in a Timoshenko beam for the pure bending case can be described by where u 1 , u 2 and u 3 are the x-, y-and z-components of the displacement vector u of a generic point with coordinates (x, y, z) on a beam cross section, v is the z-component of the displacement vector of a generic point along the x-axis and ϑ is the rotation angle (about the y-axis) of the cross section with respect to the z-axis. From Eqs. (2) and (4), it follows that the strain components are whereas the stress components are Upon denoting with ρ and h the mass density and the beam depth, respectively, the equilibrium equations of the beam for the dynamic problem are: where v (x, t) is the transverse deflection, ϑ (x, t) the rotation of a transverse section normal to the longitudinal axis, E the Young's modulus, G the shear modulus, A the cross-sectional area, I the moment of inertia, K s the shear correction factor and q (x) the distributed transverse load. The Clapeyron theorem yields to the following expression for the total energy: where M is the bending moment, V the shear action, χ the bending curvature, γ the shear deformation angle and l the element length. By substituting both the compatibility and the constitutive equations into (9), one can write the total potential energy as quadratic functional in two fields as follows: The continuous weak formulation for the Timoshenko beam element is obtained by applying the principle of stationarity of total potential energy; one can write the problem as where δv and δϑ play the role of test functions. The discrete problem takes the form At this point, the displacement v h and the rotation ϑ h are written as linear combinations of basis functions that span, respectively, the finite-dimensional spaces V 1h = span ψ (1) where ψ (1) i are the shape functions associated with the transverse displacement, whereas ψ (2) i are relative to the nodal rotations. The approximations can be written as Substitution of (14), δv = ψ (1) i (12) leads to the following compact form of the problem: where the stiffness matrix is given by and the mass matrix is where As suggested in various papers [5][6][7], the choice of Hermite cubic interpolating functions for the transverse deflection and quadratic interpolating functions for the rotation leads to the following shape functions: where the parameters Ω, ζ and η are: EI GAK s l 2 , 1 + 12Ω, x l .
Upon introducing the parameters ξ = 1 − 6Ω and κ = 1 + 3Ω, the substitution of (19) into (16) leads to while, by substituting (19) into (17), we obtain several parts that compose the element mass matrix given as follows: The element mass matrix is given by the combination of the matrices described above: It is easy to observe that in the thin beam limit, i.e., Ω → 0 (ζ → 1, ξ → 1, κ → 1), (20) reduces to the Euler-Bernoulli elementary stiffness matrix. This choice of finiteelement model based on the two-component form of the Timoshenko beam theory is powerful for the dynamic analysis; in fact, the mass matrix ensures a superconvergent method [6].

Damage Detection: Direct Problem
Consider the following abstract problem: Find x such that where x is the design variable vector, p is the set of parameters, which the solution depends on and F is the objective function. When x is the unknown and F and p are given, we have to do with a direct problem. So, in our case, the direct problem consists in the analysis of the variations in behavior, for both the static and the dynamic problems, due to the presence of a damaged beam segment. First of all, it is very important to define the suitable number of parameters needed to describe accurately the presence of diffused cracking; as suggested in [22], if damage affects a diffused zone of a beam, it can be described by only three parameters that are the position X D , the extension L D of the cracked zone and the stiffness EI D of the damaged beam, whereas the virgin one is denoted by EI U . So, the physical problem can be described as shown in Fig. 1. In computational terms, it may be appropriate to work with a dimensionless problem: Through the definition of the following nondimensional parameters it is possible to study the problems shown in Figs. 3 and 10. As one can see in Fig. 2, the beam length has been split into three segments, each one characterized by a different axial coordinate, say ξ 1 , ξ 2 and ξ 3 , that are, respectively, defined as The axial coordinate ξ i (with i = 1, 2, 3) varies in a range between 0 and 1. This subdivision permits to treat separately the three segments and to assemble the global stiffness matrix by taking into account the lower stiffness of the second segment. Furthermore, as suggested in [21,22,24,25], it is expedient to define a relationship among the eigenvalues Λ and the damage parameters, which takes the form for the physical problem, and for the non-dimensional one. Starting from the modal finite-element analysis of both the undamaged and the damaged beams, one can straightforwardly evaluate the function (25), or (26), through the computation of the ratio where superposed U and D, respectively, indicate undamaged and damaged states, Λ i is the ith system eigenvalue and ω i ≡ √ Λ i is the ith natural frequency.

Numerical Example: Simply Supported Beam
The first numerical example consists in the simply supported beam subdivided into 400 elements shown in Fig. 3, subjected to an uniform transverse load whose geometric and mechanical properties are listed in Table 1.  What one would expect from the numerical computation of the problem is an increase in the deformability of the beam due to the presence of a cracked zone. On the basis of these considerations, in the remainder of this section the results obtained for the static and the modal analysis are shown.
From Fig. 4, it is possible to observe that the maximum deflection of the damaged beam is obviously greater than those obtained with the undamaged one and that it is not located in the middle span of the beam, but it is shifted toward the cracked zone.
From Fig. 5, the increase in rotation to which the damaged beam is subjected is sharply evident: In the particular case of a concentrated damage, one can note that, in correspondence of the cracked zone, the trend of the nodal rotations undergoes a jump that can be taken into account by introducing an equivalent torsional spring. As suggested in [22], such torsional spring, characterized by a stiffness K , represents both the extension and the magnitude of the damage. As observed from results of the static problem, also for the free vibrations analysis whose results are shown in Figs. 6, 7 and 8, it is possible to note that greater variations in the response were found in terms of rotations, where the presence of a cracked zone causes an increase in the nodal rotations. Concerning the displacements field, the effects of damage are only lightly visible. In this case, we set the same finite-element size for the three segments with the aim to obtain a more clear and representative result.
In this paper, the first three frequencies are studied and the curves shown in Fig. 9 are obtained by varying both the damage position x D and the non-dimensional magnitude β.
The variation of frequencies, as expected, assumes the maximum values when the damage is located in the maximum or minimum values of the undamaged frequency modes. Furthermore, variation of frequencies increases with the damage magnitude and, obviously, the functions plotted above are equal to zero if β = 0, i.e., if the beam is not damaged.
The trends of the curves shown above are equal to those obtained in [22] for the same problem, i.e., the case of a simply supported beam via analytical solution.

Numerical Example: Cantilever Beam
A necessary condition for an inverse problem solver to work properly is of course its capability to address successfully the analysis problem (direct problem). To show the capabilities of the proposed approach, a few numerical studies concerning the cantilever in Figure 10 are then reported next. Geometric and mechanical properties of the structure are listed in Table 2, and 150 elements have been used for the discretization of the beam.
As to the damage parameters, we set a position X D = L/5 = 0.6, an extension L D = L/10 = 0.3 and a non-dimensional magnitude β = 0.5.
In this case, the effects of the presence of a cracked zone on the static analysis are visible only if the damage position is sufficiently near to the clamped end. Since a suitable damage position is defined, also for this example one can observe the increase in both the displacements and the rotations due to the stiffness reduction that affects the damaged segment.  In particular, from Fig. 11 it is possible to observe that the differences in the static response in terms of displacements begin to be visible exactly in correspondence of the cracked zone and that they grow along the beam axis.
Concerning nodal rotations, the static response for both the damaged and the undamaged beam is the same up to the start of the cracked segment, then nodal rotations of the first one suffer a rise, and finally, after the end of the damaged zone, the two solutions keep about parallel along the beam axis. If the position of the damaged zone is near to the free end of the beam, the gap between solutions is much smaller than that obtained by the proposed numerical example (Fig. 12).
Concerning free vibration analysis, also for the case of a cantilever beam the greater variations in behavior are visible from the results in terms of nodal rotations, whereas displacements are about the same for both the undamaged and the damaged beam.  In particular, nodal displacements for the first two modes (Figs. 13, 14) are approximately coincident, while in the third mode shown in Fig. 15, one can note some differences in correspondence of the cracked zone. On the contrary, the variations in behavior of the rotation field are visible in all of the three modes considered in this example. Also for this structure, in the free vibration analysis it is assumed that the finite-element length adopted for the discretization of the three segments must be the same with the aim to provide more representative results. As one can observe from Fig. 16, variations of frequencies, given from the evaluation of (27), assume their maximum value where the flexural curvature is greater and they decrease where the curvature is lower.
Furthermore, also in this numerical example, it is possible to note that variations of frequencies are strongly dependent on the damage magnitude; in fact, they grow with this last one and if β → 0, then the non-dimensional frequency variation Δω i /ω U i → 0 (i = 1, 2, 3).

Damage Detection: Inverse Problem
Considering now the problem (22), when F and x are given and p is the unknown, we have to do with an inverse problem. The damage identification procedure considered in this paper is based on the comparison between the frequencies: This method, as suggested in [21,22,24], is well known as response quantities procedure. The aim of the inverse problem is to define the damage parameters by using the following input data: -frequencies of the beam in the undamaged condition -experimental frequencies of the real beam subjected to damage.
It is important to specify that, in this paper, experimental data are given through the development of the direct problem: In such a way, one can control whether or not the damage parameters are really found and test the robustness of the method.

Numerical Optimization
This is a typical case of nonlinear constrained optimization (NCO) because the objective function is nonlinear, the damage parameters are subjected to linear inequalities and they are characterized by both lower and upper bounds. So, the general structure of a NCO problem is the following (vector inequalities are intended componentwise): In our case, the objective function used, which, as previously mentioned, is based on the comparison between experimental and analytical frequencies, is where ω D i is the experimental value of the ith frequency of the real beam subjected to damage,ω D i (x D , b D , β) is the analytical expression of the ith frequency defined as a function of the damage parameters and ω U i is the frequency value of the undamaged beam. With reference to (28), and taking into account of the non-dimensional parameters (23), the linear inequalities b D /2 ≤ x D ≤ 1 − b D /2 can be expressed in matrix form as Furthermore, the non-dimensional parameters (23) are bounded as follows: and in vectorial form assume the following form Finally, the linear equalities constraints and the nonlinear equalities/inequalities constraints in (28) are not present in this optimization problem. Since the objective function is non-convex, the traditional optimization methods (such as barrier m., penalty m. and trust region m.) allow to define the damage parameters only if the initial point x 0 is sufficiently close to the global minimum. In constrained optimization, the general aim is to transform the problem into an easier subproblem that can then be solved and used as the basis of an iterative process. To find some way, so that damage identification is successful for any choice of the initial point, a two-phase procedure has been considered in this paper. So, the objective function directly derived by Eq. 28 becomeŝ for prescribed values of both the damage extension b D and the coefficient of stiffness reduction β, and it depends only on the damage position x D . It is important to note that, obviously, the prescribed values of the two parameters b D and β must be respectful of the linear inequalities and the bounds. In such a way, the values of the damage parameters that minimize the differences between experimental and analytical frequencies are uniquely defined by using only the first three frequencies of the beam. The algorithm adopted is based on the active-set method that focuses on the solution of the Karush-Kuhn-Tucker (KKT) equations. The KKT equations are necessary conditions for optimality for a NCO problem. Denoted with n the number of damage parameters and with m the number of constraints, referring to the problem (28), the Kuhn-Tucker theorem says us that if the constraints are qualified in x * i (i = 1, . . . , n), then there exists a vector of multipliers λ * j ( j = 1, . . . , m), known as Kuhn-Tucker multipliers, such that x * i and λ * j are solutions of the following system of optimality conditions: From the m conditions (items d, e, f), it is clear that if a multiplier λ * j is positive, then it must be true that b j − A ji x * j = 0 or, in other words, that this constraint is active. It follows that in the conditions (items a, b, c) only the active constraints in x * i intervene. The procedure of damage detection described in this section has been integrated in an in-house code. The integration of the damage procedure is truly necessary to solve automatically, at each step of optimization problem, the numerical finite-element analysis.

Numerical Examples
The numerical examples considered for the optimization procedure are relative to the structures analyzed in the previous section. So, both geometric and mechanical properties of this structure are those listed in Tables 1 and 2, whereas the values of the damage parameters vary as shown in Table 3. By applying the procedure described in Sect. 4.1 to the input data listed in Table 3, we obtain the results shown in Figs A very small finite-element length has been chosen for this investigation that allowed us to approach the absolute minimum of the quadratic error function (that theoretically tends to zero when the number of elements goes to infinity).
As to Figs. 17, 18, 19 and 20 denoted by (b), which contain the interpolating surface of the objective function f (x D , b D , β) for a prescribed value of x min D that is the minimum value obtained by the optimization procedure, the global minimums are indicated with a red square.
In practical cases, it is hard to get extremely accurate frequencies and this can cause a wrong damage identification: To avoid such type of problem, one may use a greater number of frequencies to search for the absolute minimum. In spite of this, the use of higher frequencies introduces irregularities in the objective function that are likely to make damage detection more complex. In fact, by observing Figs. 17, 19 and 20 it is possible to note that the objective function presents more than one minimum and this calls for higher precision in on-site frequency measurement. Therefore, the  fundamental requirements for damage detection in the real case are both the high quality of frequencies and a great number of experimental measurements. Finally, the algorithm implemented for the development of the inverse problem is rather simple but considerably onerous in computational terms because it scans the whole structure to find damage parameters values that minimize the difference between experimental and analytical frequencies. It follows that a damage identification problem developed with such numerical optimization technique and applied to a complex structure may require a considerable execution time, even though it guarantees the definition of the global minimum.

Conclusions
After a brief introduction of the governing strong form, the Timoshenko interdependent interpolation element has been derived using a weak formulation setting. As for transverse displacement and nodal rotations, Hermite interpolation functions and quadratic functions have been used, respectively.
By using the interdependent interpolation element and with reference to a beam subjected to a single diffused cracked zone, the direct problem that consists in the study of the change in behavior of the beam due to the presence of damage is investigated. It results that the damaged beam presents a greater deformability than the undamaged one and that the differences are visible, especially in terms of nodal rotations for both the static and the free vibrations analysis. Furthermore, the non-dimensional variations of the first three frequencies of the beam are studied by varying the values of the damage parameters.
Starting from the first three frequencies obtained by the development of the direct problem, the inverse problem that consists in the definition of the unknown damage parameters of a structure is presented through the implementation of an optimization algorithm organized in a two-phase procedure by which we minimize an objective function based on the difference between experimental and analytical frequencies. This algorithm focuses on the solution of the Karush-Kuhn-Tucker equations and guarantees the definition of the global minimum in spite of its great computational cost.
Through the analysis of two numerical examples, it is demonstrated that only the first three frequencies are required to find the global minimum of the objective function for the experimental study, whereas in the real case it is necessary to use more than three frequencies because their measurements are often affected by errors that can compromise the accuracy of the solution.