Vibration Analysis of Cracked Beams Using the Finite Element Method

1. Introduction

Most of the members of engineering structures operate under loading conditions, which may cause damages or cracks in overstressed zones. The presence of cracks in a structural member, such as a beam, causes local variations in stiffness, the magnitude of which mainly depends on the location and depth of the cracks. These variations, in turn, have a significant effect on the vibrational behavior of the entire structure. To ensure the safe operation of structures, it is extremely important to know whether their members are free of cracks, and should any be present, to assess their extent. The procedures often used for detection are direct procedures such as ultrasound, X-rays, etc. However, these methods have proven to be inoperative and unsuitable in certain cases, since they require expensive and minutely detailed inspections [1]. To avoid these disadvantages, in recent decades, researchers have focused on more efficient procedures in crack detection using vibration-based methods [2]. Modelling of a crack is an important aspect of these methods.

The majority of published studies assume that the crack in a structural member always remains open during vibration [3-7]. However, this assumption may not be valid when dynamic loadings are dominant. In this case, the crack breathes (opens and closes) regularly during vibration, inducing variations in the structural stiffness. These variations cause the structure to exhibit non-linear dynamic behavior [8]. The main distinctive feature of this behavior is the presence of higher harmonic components. In particular, a beam with a breathing crack shows natural frequencies between those of a non-cracked beam and those of a faulty beam with an open crack. Therefore, in these cases, vibration-based methods should employ breathing crack models to provide accurate conclusions regarding the state of damage. Several researchers [9-11] have developed breathing crack models considering only the fully open and fully closed crack states. However, experiments have indicated that the transition between these two crack states does not occur instantaneously [12]. In reference [13] represented the interaction forces between two segments of a beam, separated by a crack, using time-varying connection matrices. These matrices were expanded in Fourier series to simulate the alternation of a crack opening and closing. However, the implementation of this study requires excessive computer time. In references [14, 15] considered a simple periodic function to model the time-varying stiffness of a beam. However, this model is limited to the fundamental mode, and thus, the equation of motion for the beam must be solved.

A realistic model of a breathing crack is difficult to create due to the lack of fundamental understanding about certain aspects of the breathing mechanism. This involves not only the identification of variables affecting the breathing crack behavior, but also issues for evaluating the structural dynamic response of the fractured material. It is also not yet entirely clear how partial closure interacts with key variables of the problem. The actual physical situation requires a model that accounts for the breathing mechanism and for the interaction between external loading and dynamic crack behavior. When crack contact occurs, the unknowns are the field singular behavior, the contact region and the distribution of contact tractions on the closed region of the crack. The latter class of unknowns does not exist in the case without crack closure. This type of complicated deformation of crack surfaces constitutes a non-linear problem that is too difficult to be treated with classical analytical procedures. Thus, a suitable numerical implementation is required when partial crack closure occurs.

In reference [16] constructed a lumped cracked beam model from the three-dimensional formulation of the general problem of elasticity with unilateral contact conditions on the crack lips. The problem of a beam with an edge crack subjected to a harmonic load was considered in [17]. The breathing crack behavior was simulated as a frictionless contact problem between the crack surfaces. Displacement constraints were applied to prevent penetration of the nodes of one crack surface into the other crack surface. In reference [18] studied the problem of a cantilever beam with an edge crack subjected to a harmonic load. The breathing crack behavior was represented via a frictionless contact model of the interacting surfaces. In [19] studied the effect of a helicoidal crack on the dynamic behavior of a rotating shaft. This study used a very accurate and simplified model that assumes linear stress and strain distributions to calculate the breathing mechanism. The determination of open and closed parts of the crack was performed through a non-linear iterative procedure.

This chapter presents the vibrational behavior of a beam with a non-propagating edge crack. To treat this problem, a two-dimensional beam finite element model is employed. The breathing crack is simulated as a full frictional contact problem between the crack surfaces, while the region around the crack is discretized into a number of conventional finite elements. This non-linear dynamic problem is solved using an incremental iterative procedure. This study is applied for the case of an impulsive loaded cantilever beam. Based on the derived time response, conclusions are extracted for the crack state (i.e. open or closed) over the time. Furthermore, the time response is analyzed by Fourier and continuous wavelet transforms to show the sensitivity of the vibrational behavior for both a transverse and slant crack of various depths and positions. Comparisons are performed with the corresponding vibrational behavior of the beam when the crack is considered as always open. To assess further the validity of this technique, the quasi-static problem of a three-dimensional rotating beam with a breathing crack is also presented. The formulation of this latter problem is similar to the former one. The main differences are: the inertia and damping terms are ignored, any possible sliding occurs in two dimensions and the iterative procedure is applied to load instead of time increment. The flexibility of the rotating beam and the crack state over time are presented for both a transverse and slant crack of various depths. The validation of the present study is demonstrated through comparisons with results available from the literature.

2. Finite element formulation

In the following, both a two and three-dimensional beam models with a non-propagating surface crack are presented. For both models the crack surfaces are assumed to be planar and smooth and the crack thickness negligible. The beam material properties are considered linear elastic and the displacements and strains are assumed to be small. The region around the crack is discretized into conventional finite elements. The breathing crack behavior is simulated as a full frictional contact problem between the crack surfaces, which is an inherently non-linear problem. Any possible sliding is assumed to obey Coulomb’s law of friction, and penetration between contacting areas is not allowed. The non-linear dynamic problem is discussed for the two-dimensional model and the corresponding quasi-static for the three-dimensional model. Both problems are solved utilizing incremental iterative procedures. For completeness reasons, the contact analysis in three dimensions and the formulation and solution of the non-linear dynamic problem are presented below.

1. Two-dimensional model

Figure 1 illustrates a two-dimensional straight cantilever beam with a rectangular cross-section $$b\times h$$ and length $$L$$ . A breathing crack of depth $$a$$ exists at position $$L_c$$ . The crack is located at the upper edge of the beam and forms an angle $$\theta$$ with respect to the $$x-$$ axis of the global coordinate system $$x,y$$ . An impulsive load is applied transversally at point $$A$$ (Figure 1).

Figure 1.

Cracked two-dimensional beam model.

In the finite element method (FEM) framework, the equilibrium equation governing the dynamic behavior of the model is:

$$M\ddot{U}+C\dot{U}+KU=R$$ E1

where $$M$$ , $$C$$ , and $$K$$ are the mass, damping, and stiffness matrices, respectively. The time-dependent vectors $$\ddot{U}$$ , $$\dot{U}$$ , $$U$$ , and $$R$$ denote the nodal accelerations, velocities, displacements, and external forces, respectively, in terms of a global Cartesian coordinate system $$x,y$$ .

1. Three-dimensional model

Figure 2 depicts a three-dimensional cantilever beam with length $$2L$$ and circular cross-section of radius $$R$$ . A breathing crack of depth $$a$$ exists at the middle of the beam. The crack has either straight or curved front (Figures 2b and 2c). The slant crack forms angles $${\theta }_y$$ , $${\theta }_z$$ with respect to $$(x,y)$$ and $$(x,z)$$ planes, respectively. Two different load cases are separately applied at the tip of the cantilever beam, i.e., twisting moment $$T$$ and bending moment $$M$$ , respectively. The bending moment is applied in several aperture angles $$M=M(\text{φ})$$ , in order to simulate a rotating load on a fixed beam. The components $$M_y=M_y(\text{φ})$$ and $$M_z=M_z(\text{φ})$$ of the bending moment $$M(\text{φ})$$ , along the directions of axes $$y$$ and $$z$$ , respectively, are functions of the aperture angle.

Figure 2.

Cracked three-dimensional beam; (a) cracked beam subjected to bending and twisting moments; (b) straight-front crack, and (c) curved front crack.

The equilibrium equation governing the quasi-static behavior of the model is:

$$KU=R$$ E2a

where $$K$$ is the stiffness matrix, $$U$$ denotes the nodal displacements and $$R$$ represent the external forces, respectively, in terms of a global Cartesian coordinate system $$x,y,z$$ .

The equilibrium equation governing the quasi-static behavior of the model is:

$$KU=R$$ E2b

where $$K$$ is the stiffness matrix, $$U$$ denotes the nodal displacements and $$R$$ represent the external forces, respectively, in terms of a global Cartesian coordinate system $$x,y,z$$ .

1. Crack modelling

Considering the three-dimensional model, the crack is composed of two surfaces, which intersect on the crack front. Parts of these two surfaces may come into contact on an interface. The size of the interface can vary during the interaction between the load and the structure, but the interface is usually comprised of two parts, i.e., an adhesive part and a slipping part, depending on the friction conditions maintained between the contacting surfaces. In the open crack state, the corresponding part of the crack surface is subjected to traction-free conditions. The so-called slave–master concept that is widely used for the implementation of contact analysis is adopted in this work for prediction of the crack-surface interference. One of the two crack surfaces is considered as the master surface, with the other as the slave. Both master and slave crack surfaces are defined by the local coordinate systems $$({}^{\text{J}}x{}_1,{}^{\text{J}}x{}_2,{}^{\text{J}}x{}_3)$$ , with $$\text{J}=\text{I}$$ for master surface and $$\text{J}=\text{II}$$ for slave surface. The axes $${}^{\text{J}}x{}_3$$ define the direction of the unit outward normal vector of the corresponding surfaces. The nodes that belong to the master and slave surfaces are called the master and slave nodes, respectively. Contact segments that span master nodes cover the contact surface of the structure. Therefore, the above problem can be regarded as contact between a slave node and a point on a master segment. This point may be located at a node, an edge, or a point of a master segment. A slave node makes contact with only one point on the master segment, but one master segment can make contact with one or more slave nodes at each time. For each contact pair, the mechanical contact conditions are expressed in a local coordinate system in the direction of the average normal to the boundaries of the bodies. Symbols $$u_i$$ and $$R_i$$ , $$i=1,2,3$$ denote nodal displacement and force components, respectively, defined on the local coordinate systems $$({}^{\text{J}}x{}_1,{}^{\text{J}}x{}_2,{}^{\text{J}}x{}_3)$$ , $$\text{j}=\text{I},\text{II}$$ . The subscripts that indicate nodal numbers are dropped for simplicity from this point forward.

Recalling the equilibrium condition, the force between the components is always expressed by the following equations:

$${}^{\text{I}}R{}_i+{}^{\text{II}}R{}_i=0,i=1,2,3$$ E3

In the open crack state, the following traction-free conditions are held between the components:

$${}^{\text{I}}R{}_i={}^{\text{II}}R{}_i=0,i=1,2,3$$ E4

From the definition of adhesion, the displacement components on the corresponding crack surfaces are interconnected by the equations:

$${}^{\text{I}}u{}_i+{}^{\text{II}}u{}_i=0,i=1,2$$ E5

When an initial gap $$g^0$$ exists in the normal direction between the master and slave nodes of the corresponding node pair, the displacement component along the normal direction is:

$${}^{\text{I}}u{}_3+{}^{\text{II}}u{}_3=g^0$$ E6

The slip state does not prohibit the existence of a gap between the crack surfaces, so equation (5) is still valid in this case. However, the tangential force component is defined in terms of friction as:

$${}^{\text{I}}R{}_i\pm \mu {}^{\text{I}}R{}_3=0,i=1,2$$ E7

where $$\mu$$ is the coefficient of Coulomb friction. Concerning the corresponding two-dimensional contact analysis, the aforementioned approach is straightforwardly used neglecting one of the three directions.

1. Incremental iterative procedure

The simulation of the breathing crack behavior as a full frictional contact problem constitutes the present study as a non-linear dynamic problem. In the FEM framework, a non-linear dynamic problem described by an equation such equation (1) is solved using an implicit direct integration scheme [20]. According to this method, the solution time interval of interest $$[0,T]$$ is subdivided into $$N$$ equal time increments $$\text{Δ}t$$ , where $$\text{Δ}t=T/N$$ . The variation of accelerations, velocities, and displacements within the time increment has a certain form and depends on the type of time integration scheme. Approximate solutions of equation (1) are sought at times $$0,\text{Δ}t,2\text{Δ}t,\dots ,t,t+\text{Δ}t,\dots ,T$$ . The calculations performed to obtain the solution at time $$t+\text{Δ}t$$ require that the solutions at previous times $$0,\text{Δ}t,2\text{Δ}t,\dots ,t$$ are known. The initial conditions of accelerations, velocities, and displacements at time zero are also required. Thus, equation (1) is evaluated at time $$t+\text{Δ}t$$ as:

$$M{}^{t+\text{Δ}t}\ddot{U}+C{}^{t+\text{Δ}t}\dot{U}+K{}^{t+\text{Δ}t}U={}^{t+\text{Δ}t}R$$ E8

where the left-hand subscripts denote the time.

The solution to this non-linear problem requires an iterative procedure. Employing the modified Newton-Raphson iteration method [20], the displacement vector $${}^{t+\text{Δ}t}U{}^{(k)}$$ at time $$t+\text{Δ}t$$ and iteration $$k$$ is given by:

$${}^{t+\text{Δ}t}U{}^{(k)}={}^{t+\text{Δ}t}U{}^{(k-1)}+\text{Δ}{U}^{(k)}$$ E9

while equations (8) are written as:

$$M{}^{t+\text{Δ}t}\ddot{U}{}^{(k)}+C{}^{t+\text{Δ}t}\dot{U}{}^{(k)}+{}^{t}K{}_T\text{Δ}{U}^{(k)}={}^{t+\text{Δ}t}R-{}^{t+\text{Δ}t}F{}^{(k-1)}$$ E10

where the right-hand subscripts in brackets represent the iteration number, with $$k=1,2,3,\dots$$ . The symbols $${}^{t}K{}_T$$ , $${}^{t+\text{Δ}t}F{}^{(k-1)}$$ and $$\text{Δ}{U}^{(k)}$$ denote the tangent stiffness matrix, the nodal force vector, which is equivalent to the element stresses, and the incremental nodal displacement vector, respectively. The iterative method is called the modified Newton–Raphson method, since the tangent stiffness matrix is not calculated in every iteration, which is the case for the full method [20].

Employing for example an implicit time integration scheme, formulas are implemented that relate the nodal accelerations, velocities, and displacement vectors at time $$t+\text{Δ}t$$ to those at previous times. Considering these formulas, equation (10) can be written in the following form [20]:

$${}^t\stackrel{⌢}{K}{}_T\text{Δ}{U}^{(k)}=\text{Δ}{\stackrel{⌢}{R}}^{(k-1)}$$ E11

where the matrix $${}^t\stackrel{⌢}{K}{}_T$$ is a function of the tangent stiffness matrix, mass matrix, and damping matrix, while the vector $$\text{Δ}{\stackrel{⌢}{R}}^{(k-1)}$$ contains the nodal force vector and contributions from the inertia and damping of the system. In the first iteration, the vectors $$\text{Δ}{\stackrel{⌢}{R}}^{(0)}$$ and $${}^{t+\text{Δ}t}U{}^{(0)}$$ are equal to the corresponding vectors of the last iteration at the previous time. In each iteration, the latest estimates of displacements are used to evaluate the vector $$\text{Δ}{\stackrel{⌢}{R}}^{(k-1)}$$ . Then, the incremental displacements $$\text{Δ}{U}^{(k)}$$ are obtained by solving equations (11), while the nodal displacements $${}^{t+\text{Δ}t}U{}^{(k)}$$ are derived from equations (9). The iteration proceeds until the nodal displacements vector of the last iteration $${}^{t+\text{Δ}t}U{}^{(k)}$$ are approximately equal to the corresponding vector of the previous iteration $${}^{t+\text{Δ}t}U{}^{(k-1)}$$ . This convergence criterion is expressed as

$$\left|\frac{{}^{t+\text{Δ}t}U{}^{(k)}-{}^{t+\text{Δ}t}U{}^{(k-1)}}{{}^{t+\text{Δ}t}U{}^{(k)}}\right|\le \epsilon$$ E12

where $$\epsilon$$ is a small numerical quantity.

Considering that the problem has been solved for time $$t$$ , and consequently, vectors $${}^{t}U$$ and $${}^{t}R$$ are known for the entire structure. To determine the corresponding displacement and forces vectors at time $$t+\text{Δ}t$$ , the equations (3)-(7) are written in incremental form as following:

$${}^{t+\text{Δ}t}({}^{\text{I}}\text{Δ}{R}_i)+{}^{t+\text{Δ}t}({}^{\text{II}}\text{Δ}{R}_i)=0,i=1,2$$ E13

$${}^{t+\text{Δ}t}({}^{\text{I}}\text{Δ}{R}_i)=-{}^t({}^{\text{II}}R{}_i),i=1,2$$ E14

$${}^t({}^{\text{I}}u{}_1)+{}^{t+\text{Δ}t}\text{(}{}^{\text{I}}\text{Δ}{u}_1)={}^t({}^{\text{II}}u{}_1)+{}^{t+\text{Δ}t}\text{(}{}^{\text{II}}\text{Δ}{u}_1)$$ E15

$${}^t({}^{\text{I}}u{}_2)+{}^{t+\text{Δ}t}\text{(}{}^{\text{I}}\text{Δ}{u}_2)={}^t({}^{\text{II}}u{}_2)+{}^{t+\text{Δ}t}\text{(}{}^{\text{II}}\text{Δ}{u}_2)-g^0$$ E16

$${}^t({}^{\text{I}}R{}_1)+{}^{t+\text{Δ}t}\text{(}{}^{\text{I}}\text{Δ}{R}_1)=\pm \mu ({}^t({}^{\text{I}}R{}_2)+{}^{t+\text{Δ}t}\text{(}{}^{\text{I}}\text{Δ}{R}_2))$$ E17

Assumption	Decision
Open	Contact
Open	$${}^{\text{II}}\Delta u_3{}^m-{}^{\text{I}}\Delta u_3{}^m{}^{\text{I}}u{{}_3}^{m-1}-{}^{\text{II}}u{{}_3}^{m-1}+g^0$$	$${}^{\text{II}}\Delta u_3{}^m-{}^{\text{I}}\Delta u_3{}^m\le {}^{\text{I}}u{{}_3}^{m-1}-{}^{\text{II}}u{{}_3}^{m-1}+g^0$$
Contact	$${}^{\text{I}}f{{}_3}^{m-1}+{}^{\text{I}}\Delta f_3{}^m\ge 0$$	$${}^{\text{I}}f{{}_3}^{m-1}+{}^{\text{I}}\Delta f_3{}^{m}0$$
	Adhesion	Slip
Adhesion	$$\left|{}^{\text{I}}f{{}_i}^{m-1}+{}^{\text{I}}\Delta f_i{}^m\right|\left|\mu ({}^{\text{I}}f{{}_3}^{m-1}+{}^{\text{I}}\Delta f_3{}^m)\right|$$, $$i=1,2$$	$$\left|{}^{\text{I}}f{{}_i}^{m-1}+{}^{\text{I}}\Delta f_i{}^m\right|\ge \left|\mu ({}^{\text{I}}f{{}_3}^{m-1}+{}^{\text{I}}\Delta f_3{}^m)\right|$$, $$i=1,2$$
Slip	$$({}^{\text{I}}f{{}_i}^{m-1}+{}^{\text{I}}\Delta f_i{}^m)({}^{\text{I}}\Delta f_i{}^m-{}^{\text{II}}\Delta f_i{}^m)0$$, $$i=1,2$$	$$({}^{\text{I}}f{{}_i}^{m-1}+{}^{\text{I}}\Delta f_i{}^m)({}^{\text{I}}\Delta f_i{}^m-{}^{\text{II}}\Delta f_i{}^m)\le 0$$, $$i=1,2$$

Table 1.

Definition of contact status.

For reasons of simplicity, the iteration number has been omitted from equations (13)-(17). However, the formulation given below is repeated for all iterations. These equations are transformed to the global Cartesian coordinate system $$x,y,z$$ and are then embedded and rearranged into equation (11). To determine the corresponding nodal displacements at time $$t+\text{Δ}t$$ , the contact conditions must first be satisfied. Therefore, the iterative procedure employed must be applied by initially using the convergent contact status (union of the adhesive, slipping and open parts of the crack surface) of the previous time $$t$$ . The procedure initially assumes that the coplanar and normal incremental force components for a master surface at time $$t+\text{Δ}t$$ are zero. Accurate values of the incremental forces can be estimated via the iterative procedure. The contact state for every node pair is examined according to Table 1. This table describes criteria to check whether violations involving geometrical compatibility and force continuity have occurred. Where necessary, appropriate changes from open to contact or from adhesion to slip states and vice versa are made to identify the equilibrium state of the contact conditions. The new contact condition is applied to the node pair closest to the change. If the change is from the open to the contact state, then the adhesion condition is adjusted. When the iterative procedure converges, the incremental nodal values $${}^{t+\text{Δ}t}\text{Δ}U$$ and $${}^{t+\text{Δ}t}\text{Δ}R$$ are known for the entire structure. After calculating the total nodal values, the procedure goes to the next step of the time increment and continues until the final time increment is reached. The problem solution is then attained.

3. Local flexibilities in cracked beams

A crack introduces local flexibilities in the stiffness of the structure due to strain energy concentration. Although local flexibilities representing the fracture in stationary structures are constant for open cracks, the breathing mechanism causes their time dependence. In a fixed direction, local flexibilities of rotating beams change also with time due to the breathing mechanism. Evidently, the vibrational response of a rotating cracked beam depends on the crack opening and closing pattern in one cycle.

Since the torsional and bending vibrations are dominant in rotating beams, in this chapter it is assumed that the corresponding local flexibilities are also dominant in the local flexibility matrix, neglecting the cross-coupling terms. Numerical results showed that the off-diagonal coefficients of this matrix are at least two orders of magnitude lower than the diagonal ones, and thus are considered negligible. Therefore, the presence of the crack can equivalently represented by a diagonal local flexibility matrix, independent of the crack contact conditions.

The exact relationship between the fracture characteristics and the induced local flexibilities is difficult to be determined by the strain energy approach, because, stress intensity factor expressions for this complex geometry are not available. For the computation of the local crack compliance, a finite element method was used. According to the point load displacement method, if at some node preferable lying on the tip of the beam is applied the external load vector $$\left\{Q\right\}=\left\{\begin{array}{ccc}T& M_y& M_z\end{array}\right\}{}^{{\rm T}}$$ , then at the same node the resulting rotations $$\left\{\text{Θ}\right\}=\left\{\begin{array}{ccc}{\text{Θ}}_{\text{x}}& {\text{Θ}}_y& {\Theta }_z\end{array}\right\}{}^{{\rm T}}$$ are:

$$\left[\begin{array}{ccc}{c}_x& 0& 0\\ 0& c_y& 0\\ 0& 0& c_z\end{array}\right]\left\{\begin{array}{c}T\\ M_y\\ M_z\end{array}\right\}=\left\{\begin{array}{c}{\text{Θ}}_{\text{x}}\\ {\text{Θ}}_{\text{y}}\\ {\text{Θ}}_{\text{z}}\end{array}\right\}$$ E18

In equation (18), $$c_r$$ , $$r=x,y,z$$ are the diagonal coefficients of the local flexibility matrix. Under the application of a particular load component $$Q_r$$ at the $$r$$ - direction, the above equation is then simplified as follows:

$${\Theta }_r=c_r{Q}_r$$ E19

where $${\Theta }_r$$ , $$r=x,y,z$$ is the induced rotation, and $$c_r$$ is the local flexibility component that corresponds to the particular loading mode $$Q_r$$ . Equation (19) can be used for the computation of the local flexibility coefficients, as described in the following. When the original non-cracked beam is uploaded until the load $$Q_{ro}$$ , a rotation is imposed in $$r$$ - direction, such that equation (19) gives:

$$q_{ro}=c_{ro}{Q}_{ro}$$ E20

where $$c_{ro}$$ is the flexibility of the original structure. The deformation of the cracked beam when loaded at the same node gives:

$${\Theta }_{rc}=c_{rc}{Q}_{rc}$$ E21

where $$c_{rc}$$ is the flexibility of the cracked beam, and $$Q_{rc}$$ , $${\Theta }_{rc}$$ the applied load and the resulting rotation, respectively. Between the flexibilities of the original and the fractured structure holds the condition

$$c_{rc}=c_r+c_{ro}$$ E22

where $$c_r$$ is the local flexibility due to the crack itself. Assuming that the applied load levels are of the same magnitude, i.e. $$Q_{rc}=Q_{ro}=Q_r$$ , after some manipulation, equations (19)-(23) yield the local flexibility coefficient in the $$r$$ - direction

$$c_r=\frac{{\Theta }_{rc}-{\Theta }_{ro}}{Q_r}$$ E23

Equation (23) is used to compute the coefficients of the local flexibility matrix $$[c]$$ utilizing the FEM results. Tip loads $$Q_r$$ , $$r=x,y,z$$ are applied independently, and the resulting rotations $${\Theta }_r$$ are evaluated. The rotations of the original structure $$Q_{ro}$$ , are evaluated for FEM models that do not present crack but have similar meshing with the cracked models. When fractured models are examined, FEM results are computed for several values of crack depth and different loading conditions.

4. Response analysis

The response derived from equations (9) and (11) cannot be examined directly to distinguish the breathing crack effects. For this reason, fast Fourier and continuous wavelet transforms are employed. These two popular transforms in signal analysis are briefly discussed below for reasons of completeness, and more information can be found in references [21, 22].

The fast Fourier transform (FFT) is a perfect tool for finding the frequency components in a signal of stationary nature. Unfortunately, FFT cannot show the time point at which a particular frequency component occurs. Therefore, FFT is not a suitable tool for a non-stationary signal, such as the impulsive response of the cracked cantilever beam considered in this study, which requires time-frequency representation. To overcome this FFT deficiency, the short time Fourier transform (STFT) could be adopted, which maps a signal into a two-dimensional function of time and frequency. This windowing technique analyzes only a small section of the signal at a time. However, the information about time and frequency that is obtained has a limited precision that is determined by the size of the window, which is the same for all frequencies.

Wavelet transforms are a novel and precise way to analyze signals and can overcome the problems that other signal transforms exhibit. The most important advantage of wavelet transformations is that they have changeable window dimensions. For low frequencies, the window is wide, while for high frequencies, it is narrow. Thus, maximum time frequency resolution is provided for all frequency intervals.

The continuous wavelet transform (CWT), as employed in this study, is defined mathematically as:

$$W{f}_{s,u}=\frac{1}{\sqrt{s}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}f(t)}{\psi }^{\ast }(\frac{t-u}{s})dt$$ E24

where $$f(t)$$ is the signal for analysis, $${\psi }^{\ast }(t)$$ is the complex conjugate of the mother wavelet $$\psi (t)$$ and $$s$$ and $$u$$ are real-valued parameters used to characterize the dilation and translation features of the wavelet.

The CWT has an inverse that permits recovery of the signal from its coefficient $$W{f}_{s,u}$$ and is defined as:

$$f(t)=\frac{1}{C_{\psi }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}W{f}_{s,u}\psi (\frac{t-u}{s})\frac{1}{s^2}}}dsdu$$ E25

5. Results and discussions

1. Accuracy study

1. Two-dimensional model

To demonstrate the accuracy of the presented study, a two-dimensional beam is considered with length $$L=\text{1}\text{.5m}$$ , cross-section $$a/R=1.0$$ , modulus of elasticity $$E=2.06\times {10}^{11}\text{Pa}$$ , mass density $$\rho =7650{\text{kg/m}}^{\text{3}}$$ and Poisson’s ratio $$\nu =0.29$$ . It is assumed that the beam contains a breathing crack of $$\theta ={90}^{\circ }$$ , $$L_c/L=0.5$$ and $$a/h=0.5$$ . A transverse impulse loading is applied at point $$A$$ of the beam (Figure 1). Based on the FFT of the transverse acceleration at point $$L_r/L=1$$ (Figure 1), the three lower dimensionless natural frequencies $$f_{ic}/f_{in}$$ of the cantilever beam are evaluated and quoted at Table 2. Subscript $$i=1,2,3$$ denotes the order of natural frequency, while the subscript $$j$$ represents the crack state of the beam. In particular, for $$j=c$$ the beam is cracked, while for $$j=n$$ the beam is non-cracked. It is derived from Table 2, that the results of the present study are generally close to the results of references [23] and [24]. It is noteworthy that the study of reference [23] has applied for a Timoshenko beam. Furthermore, the results from the work in [24] correspond to an internal crack of the same severity with the studied crack case. The internal crack constitutes a good approximation of a fully closed crack.

	$$f_{1c}/f_{1n}$$	$$f_{2c}/f_{2n}$$	$$f_{3c}/f_{3n}$$
Nandwana and Maiti [23]	1.000	0.992	1.014
Present study	0.986	0.992	0.998
% difference of natural frequency	1.41	-0.04	1.60
Kisa and Brandon [24]	0.980	0.925	0.999
Present study	0.986	0.992	0.998
% difference of natural frequency	-0.58	-6.71	0.08

Table 2.

Three lower natural frequencies of the cantilever beam with a crack of $$\theta =90\text{°},L_c/L=0.5$$ and $$a/h=0.5$$ .

1. Three-dimensional model

Figure 3.

Flexibility of a transverse straight front crack under subjected to traction-free conditions; (a) bending flexibility coefficient $$C_z$$ , and (b) torsional flexibility coefficient $$C_x$$ .

To demonstrate further the accuracy of the presented study, a three-dimensional beam with length $$L=\text{1}\text{.0m}$$ , radius $$R=\text{0}\text{.05m}$$ , modulus of elasticity $$E=2.06\times {10}^{11}\text{Pa}$$ , mass density $$\rho =7650{\text{kg/m}}^{\text{3}}$$ and Poisson’s ratio $$\nu =0.29$$ is considered. The beam has a gaping straight front crack of various depths (Figure 3). The beam undergoes either pure bending $$M_z$$ , i.e. $$\varphi ={90}^o$$ or twisting moment $$T$$ . For both loading cases, the dimensionless local flexibility coefficients, $$C_{r=}(E{R}^3/(1-{\nu }^2))c_r$$ with $$r=x,y,z$$ , are calculated for various values of dimensionless crack depth $$a/h$$ . The application of bending moment for this transversely fractured structure imposes the bending mode local flexibility $$C_z$$ (Figure 3a), while the twisting moment imposes only the twisting coefficient is $$C_x$$ (Figure 3b). As shown from Figures 3, the presented results are in good agreement with semi-analytical and experimental ones [25-27].

Figure 4.

Dimensionless bending compliance variation of a transversely cracked beam as function of the aperture angle; (a) straight front crack, and (b) curved front crack.

Figure 4 presents the variation of the flexibility coefficient $$C_y$$ when the transversely cracked beam is subjected to bending moment $$M=M(\varphi )$$ . The numerical results are plotted for a period that corresponds to one revolution of the beam. For better reasons of understanding, the local flexibility is plotted with phase lag $$\pi /2$$ that corresponds to the application only of the $$M_z$$ component of bending moment or $$\text{φ}={90}^o$$ . That is, the evolution of the local flexibility is progressing from the full open condition of the crack. The same figure illustrates theoretical [28], experimental [29] and numerical [30] results from the literature for straight-fronted edge crack and crack depth $$a/R=1.0$$ . As it is shown, the values of local flexibility for crack depth $$a/R=1.0$$ , from FEM analysis are smaller than the other methods, when the crack is open. The experimental results are obtained for notched beam which is more flexible than the cracked beam. As the crack closes the FEM results are greater than the others. This is explained by the crack contact area (Figure 5), which differs and is generally smaller than the area assumed by the theoretical approaches for the same edge orientation. The impossibility of these approximations to predict the crack closure correctly is the main reason that these yield comparable results with the present ones only on the regions of partially opening portion of the crack, i.e., between $$L_c/L=0.5$$ or $$a/h=0.75$$ . The variation of the coefficient $$C_y$$ depends significantly on the crack depth. For small crack depths, the crack does not open regularly once per revolution, but contact is observed twice per revolution. The second contact state yields smaller compliance than the regular contact.

1. Crack contact state

Figure 5 illustrates predictions of the crack closure portion for a transverse crack with depth $$a/R=0.8$$ , when the three-dimensional beam is loaded in bending moment only, versus the aperture angle. Figure 5a shows the crack closure evolution for a straight front crack, and Figure 5b for a curved front crack, respectively. As the aperture angle increases, the crack closure portion increases in both cases examined. On full load reversal, a very small portion of the crack surface along the crack front remains always open. The shape of the contact surface and its portion clearly depends on the shape of the crack. This fact is expected to impose differences between the local flexibilities of different crack shapes. When the crack is slant, smaller portions of the crack surfaces are in contact.

Figure 5.

Evolution of contact area between the crack surfaces of a transverse crack under bending loading; (a) straight front crack, and (b) curved front crack.

Figure 6 depicts this situation for cracks having slope $${\theta }_z$$ with respect to the z-axis ( $${\theta }_y=0^o$$ ), loaded in bending $$M_z$$ (or $$\text{φ}={90}^o$$ ). The crack slope $${\theta }_z$$ seems to affect slightly the size and shape of the contact area between the crack surfaces for both of the crack fronts under investigation. The cases of crack orientation subjected to bending load that are not illustrated here present similar behavior with respect to the crack surface contact. For cracks subjected to twisting moment, the contact area does not generally change during the shaft rotation. Figure 7 depicts this fact when the crack slope is defined by the angles $${\theta }_y={45}^o$$ and $${\theta }_z=0^o$$ .

Figure 6.

Contact area between crack surfaces for different crack slopes when the shaft is loaded in bending ( $${\theta }_y=0$$ , $$\text{φ}={90}^o$$ ); (a) straight front crack, and (b) curved front crack.

Figure 7.

Contact area between crack surfaces when the shaft is loaded in torsion; (a) straight-fronted crack, and (b) penny shaped crack.

Further knowledge regarding the crack contact state can be also derived studying a strip under reversing loadings. Figure 8a illustrates the distinction between the open and partially closed crack configurations as defined by the crack opening displacements (CODs). In this figure the normalized displacements of the crack surface over the strip height h, close to the crack tip is shown. It is apparent that direct bending opens the crack, although load reversing causes partial crack closure. Solid lines represent crack opening mode and dashed lines partial crack closure, respectively. Negative values in this figure represent sign convention between the axes and not penetrating crack surfaces. Numerical results show that open crack displacements are two orders of magnitude higher than those corresponding to partially closed cracks. Figure 8b shows the transverse to crack line CODs, for two cracks with θ=135° and depths $$a/h=0.8$$ and $$a/h=1$$ , respectively. Under load reversing conditions and in the region close to crack tip, opening displacements develop, as in previous case. Crack orientation affects the deformation mode, yielding smaller displacements than those developed in transverse cracks. In the case of slant cracks, tangential stresses are shown to play a significant role in the crack surface interference process, producing small opening ligaments and destroying symmetry of displacements. CODs observed in slant cracks under load reversing are two orders of magnitude smaller than the corresponding displacements observed on normal loading inducing stress intensification.

Figure 8.

CODs for (a) transverse cracks, and (b) slant cracks.

1. Time Response

Figure 9.

Time response when a crack is present at $$L_c/L=0.25$$ with depth $$a/h=0.75$$ and angle $$\theta ={90}^{\circ }$$ .

Figure 9 depicts the impulsive displacement response in the transverse direction versus time $$t$$ for the two-dimensional beam with a breathing crack of $$\theta ={90}^{\circ }$$ , $$L_c/L=0.25$$ and $$a/h=0.75$$ . Consider a time period of response that spans from $$t=0\sec$$ to $$t=0.05\sec$$ . From the deformed mesh, it is assumed that during approximately the first half of the time period (from $$t=0\sec$$ to $$t=0.02\sec$$ ) where the response sign is positive, the crack opens and closes continuously. On the contrary, during the second half of the period (from $$t=0.02\sec$$ to $$t=0.05\sec$$ ) where the sign of the response is negative, the crack is open. For this period, the response is composed of a sinusoidal waveform that corresponds to that of the non-cracked beam, and several disturbances that correspond to lower frequencies than those of the non-cracked beam. For the first part of the indicated period (where the displacement is positive), these disturbances correspond to the existence of partial closure and contact along the crack surfaces, while the second part of the period (where the displacement is negative) corresponds to crack opening deformations. The same conclusions apply for the part of the response to the right of the considered time period. As it is deduced from Figure 9, the breathing of crack causes the beam to exhibit a highly nonlinear vibrational response. As a result, the portions of the beam on the left and right of the crack vibrate with different ways at the same time. In an attempt to clarify the intricacies involved in the response, Figure 10 illustrates the deformed shape and the detail of the corresponding crack opening deformation for a beam with either a small ( $$a/h=0.25$$ ) or a deep ( $$a/h=0.75$$ ) crack of angle $$\theta ={45}^{\circ }$$ at position $$L_c/L=0.25$$ for time $$t=0.013\sec$$ at which the crack is instantaneously open. Apart from the instantaneous deformed shapes of the fractured beams, Figure 10a shows also the undeformed shape of the non-cracked beam (dotted lines) for comparison. The displacements in Figure 10a are magnified with a scale $${10}^4$$ , while in Figures 10b-10c are magnified with a scale $$6\times {10}^4$$ . It seems from Figure 10a that the instantaneous deformed shapes for the small and the deep crack approach mainly the fundamental and the second mode shapes of the non-cracked beam, respectively.

Figure 10.

(a) Instantaneous deformed shape of a beam with an inclined crack ( $$\theta ={45}^{\circ }$$ , $$L_c/L=0.25$$ ) at time $$t=0.013\sec$$ and detail of the corresponding crack opening deformation for (b) a small ( $$a/h=0.25$$ ) and (c) a deep ( $$a/h=0.75$$ ) crack.

It is also visible from Figure 10a that for the deep crack the beam portion located on the left of the crack is nearly undeformed due to the dominance of the support conditions in the vibrational behavior of the beam at the considered time. Based on the Figures 10b-10c, it is observed that the distance between crack surfaces increases as the crack depth increase. This is reasonable since the small crack reduces slightly the local stiffness at crack position. This slight change of local stiffness matrix is time dependent and varies from zero to a maximum value depending on the breathing mechanism. On the contrary, for the deep crack the range of local stiffness change is wider than that of the small crack cases. For this reason, the effect of the crack is more important. The deformed shape of beam and the corresponding detail of the crack opening deformation for the two considered crack cases are also presented in Figure 11, which is magnified as Figure 10, for time $$t=0.023\sec$$ at which the crack is instantaneously closed. It seems from Figure11a that the instantaneous deformed shape for the small crack approaches mainly the fundamental mode of the non-cracked beam. In contrast, the instantaneous deformed shape for the deep crack approaches mainly the second mode shape of the non-cracked beam. It is noticeable from Figures 11b-11c that the crack surfaces are not completely closed due to the presence of the friction, which causes the development of shear stresses along crack surfaces. This stress component causes deformation of crack surfaces and develops a mechanism of energy absorption.

Figure 11.

(a) Instantaneous deformed shape of a beam with an inclined crack ( $$\theta ={45}^{\circ }$$ , $$L_c/L=0.25$$ ) at time $$t=0.013\sec$$ and detail of the corresponding crack opening deformation for (b) a small ( $$a/h=0.25$$ ) and (c) a deep ( $$a/h=0.75$$ ) crack.

1. Frequency Response

Figure 12 depicts the FFT of the transverse acceleration response for the two-dimensional beam with a breathing crack of $$\theta ={30}^{\circ }$$ , $$L_c/L=0.25$$ , and $$a/h=0.75$$ . The corresponding FFT of transverse acceleration response with the crack always open is also plotted in this figure for comparison in frequency domain. For reasons of clarity, in this and the other figure of this subchapter, the corresponding response for the non-cracked beam is not plotted in this figure. The four vertical dash-dot lines represent the loci of the first four natural bending frequencies of the non-cracked beam. These natural frequencies are evaluated from the FFT of the transverse acceleration response and are $$0.028,0.174,0.472$$ , and the $$0.873\text{kHz}$$ . The dash-dot-dot vertical lines in the figures represent the third natural bending frequency of the beam portion located on the right of the crack. For the present crack case, this portion of the beam has a length of $$0.75L$$ , and its third natural bending frequency is $$0.810\text{kHz}$$ . This latter frequency is evaluated from the FFT of the transverse acceleration response. Figure 12 shows that the FFT of the open crack model exhibits five peaks at $$0.026,0.172,0.448,0.784$$ , and $$0.850\text{kHz}$$ . The fourth peak corresponds to the third natural bending frequency of the beam portion with length $$0.75L$$ , while the remaining peaks correspond to the first four natural bending frequencies of the non-cracked beam. Corresponding peaks appear for the breathing crack model (Figure 12). The two FFTs illustrated in Figure 12 show that the first two frequencies of the open and breathing crack models are very close and are lower than those corresponding to the non-cracked beam. Shifts are observed for the remaining three frequencies. In particular, the third, fourth, and fifth frequencies of the breathing crack model are between the frequencies of the open cracked model and the non-cracked beam, as expected [8]. The absolute percentage differences for these three frequencies are $$2.40\%$$ , $$2.61\%$$ , and $$0.93\%$$ , respectively.

Figure 12.

Frequency response of the transverse acceleration when a crack of depth $$a/h=0.25$$ and angle $$\theta ={90}^{\circ }$$ is present at $$L_c/L=0.5$$ .

Figure 13 shows the sensitivity of beam vibrational behavior in terms of crack position. It seems that the crack position affects all but the first natural frequency. The absolute percentage differences are $$12.87\%$$ , $$8.53\%$$ , $$7.29\%$$ , and $$1.47\%$$ for the second through fifth natural frequencies, respectively.

Figure 13.

Frequency response of the transverse acceleration when a breathing crack of angle $$\theta ={45}^{\circ }$$ and depth $$a/h=0.75$$ is present.

1. Time-Frequency Response

Figure 14a shows a contour map of the CWT for the transverse displacement response of the two-dimensional beam with a breathing crack of $$\theta ={90}^{\circ }$$ , $$L_c/L=0.25$$ , and $$a/h=0.75$$ . In this and following contour map, the color at each point of the $$s-u$$ domain represents the magnitude of the wavelet coefficients $$W{f}_{s,u}$$ . A lighter color corresponds to larger coefficients, and a darker color to smaller ones. The symbol $$s$$ denotes the scale space that is inversely proportional to the frequency domain, and $$u$$ stands for the time space. The contour map of Figure 14a consists of three horizontal regions. The first region extends up to $$s\approx 10$$ , the second one up to $$s\approx 40$$ , and the third occupies the rest of the map. The third region is unchanged for $$s128$$ . Each of these regions consists of a number of consecutive approximately equal vertical ridges. The number of ridges differs according to region.

Higher numbers of ridges appear in the first region while there are fewer in the third. A comparison of these three regions shows differences in the magnitude of the wavelet coefficients. Figure 14b shows the contour map of the CWT for the transverse displacement response for a fully open crack. This contour map consists of three visible horizontal regions. The first two regions are similar to those of the breathing crack model (Figure 14a) in respect to their extent and the number and extent of ridges. On the contrary, the magnitude of the wavelet coefficients appears differences. The third region consists of fewer and wider ridges than that of the breathing model (Figure 14a).

Figure 14.

CWT of the transverse displacement history when (a) a breathing crack, (b) an open crack is present ( $$a/h=0.75$$ , $$\theta ={90}^{\circ }$$ , and $$L_c/L=0.25$$ ).

Figure 15.

CWT of the transverse displacement history when a breathing crack of angle $$\theta ={45}^{\circ }$$ is present at $$L_c/L=0.5$$ ; (a) $$a/h=0.25$$ , and (b) $$a/h=0.75$$ .

The higher number of ridges in the breathing model most likely represents changes to stiffness matrix due to the breathing of the crack. Similar conclusions regarding the magnitude of the wavelet coefficients, the extent of the regions, and the number and shape of ridges per region can be made from Figure 15, which shows the sensitivity analysis with respect to crack depth. Although, the effect of the crack depth in the morphology of the Figure 15 is visible, a qualitative analysis is required to correlate the natural frequencies obtained from the FFT to the extent of the horizontal regions and the number and the extent of ridges per region. Furthermore, the magnitude of the wavelet coefficients in the domain $$s-u$$ can be extracted. This analysis correlates the natural beam frequencies obtained by the FFT to the scales of contour maps. The magnitude of the wavelet coefficients at these scales can then be found from the contour maps for the entire time solution of interest, and conclusions for crack detection techniques can be extracted.

6. Conclusions

In this chapter, finite element procedures able to approach the vibrations of a beam with a breathing crack were presented. Although the developed models were discretized into a number of conventional finite elements, the breathing was treated as a full frictional contact problem between the crack surfaces. Quasi-static and non-linear dynamic analyses were performed aiming the prediction of vibration characteristics of cracked beams. The solutions were obtained using incremental iterative procedures. The results show good agreement with experimental or theoretical ones found in the open literature. The assesment of crack contact state in the different phases of the breathing mechanism gives comprehesive answers for the local flexibility variations and concequently for the vibrational response of a cracked beam. The derived time response of the two-dimensional beam was analyzed by various integral transforms including FFT and CWT. This model was assessed for the case of a cantilever beam subjected to an impulse loading. The results show the sensitivity of vibrational behavior with respect to crack characteristics. Although a qualitative analysis was required to interpret the results exactly, this study proved that FFT and CWT can be used as supplementary tools for crack detection techniques.