Nonlinear Vibration Analysis of a Beam with a Breathing Crack

Abstract

The phenomena of sub- and super-harmonic responses make up one of the prominent nonlinear characteristics of a beam with a breathing crack. In order to fully understand the behaviors of sub- and super-harmonic resonances, it is necessary to analyze the nonlinear vibration of a beam-like structure with a breathing crack. In this study, a new stiffness model that considers the influence of the partial crack closure is proposed to model the stiffness variation of the cracked beam. Based on the finite element model of a beam with a breathing crack, the multiple-scale method is proposed to analyze the nonlinear vibration of a cracked beam subjected to harmonic excitation, and the relation between the nonlinear vibration of the cracked beam and the system parameters is obtained. An experiment is conducted to validate the analytical results. The study shows that the nonlinear responses of a beam with a breathing crack are affected by both the structural parameters and the crack parameters.

1. Introduction

Beams are one of the basic structural members and are widely used in various fields, such as aerospace, machinery, automobile and structure industries [1,2,3,4]. Most beam-like structures are subjected to dynamic loading. Due to various causes such as corrosion and temperature fluctuation, cracks may result in failure [5]. If these cracks are not detected in time, a sudden structural failure can lead to catastrophic damages [6,7]. Due to the various types of crack in different structures, currently, industries have a number of different damage-detection methods, such as ultrasonic methods [8,9,10], X-ray methods and vibration-based damage identification methods. Owing to the advancement of reliable fault detection and diagnosis techniques, vibration-based damage identification methods have been developed for various engineering applications, such as the rotating blades of engines [11], submarine pipes [12], rotating rotors [13] and atomic force microscope (AFM) probes [14].

In the early stage, most researchers assumed that cracks keep opening during vibration, and thus a linear system was adopted. In this case, the dynamic characteristics of the linear system, such as natural frequencies and mode shapes, are used to determine the crack severity. However, these linear response features are not valid when detecting fatigue cracks due to the crack breathing effect in reality [15]. For example, Gudmunson [16] experimentally investigated the crack closure effect on the vibrations of a cracked beam. The results showed that the natural frequencies measured from the experiment are different from those obtained from analysis that does not consider the effect of crack closing. Therefore, in the case of the closing crack (also known as breathing crack), the dynamic system is nonlinear because of the time-varying stiffness caused by opening and closing effects of the crack [17].

Over the last decades, many scholars have exploited the different nonlinear characteristics of beams with breathing cracks by using analytical methods [18,19,20], numerical methods [21,22] and experimental approaches [23,24]. Their studies showed that nonlinear vibration characteristics have a higher sensitivity to breathing cracks, especially the presence of sub- and super-harmonic resonances [25]. For example, Matveev et al. [18] presented an approximate analytical method to calculate the relative amplitude of the dominant harmonics of a rod-like structure with a closing crack at sub- and super-harmonic resonances, and they showed that the acceleration amplitude of the second harmonic is more sensitive than that of the fundamental harmonic. Bovsunovskii et al. [19] developed the dynamic model of a beam with a closing crack which considered the damping change due to crack propagation. The relationship of the amplitude of dominant harmonics with the force application point and sensor location at sub-resonance and super-resonance was investigated. As the crack depth increases, the super-harmonic resonance in the frequency response is distinct, meaning that the second order becomes dominant. Andreaus et al. [21,22] considered the closing crack interfaces as a contact problem and studied the forced harmonic response of the cracked beam by using the finite analysis software ADINA. The simulation results demonstrated that the phenomena of sub- and super-harmonic resonances exist in the cracked beam. Later, the same authors [24] also conducted an experimental investigation into the forced harmonic response of a cracked beam with a breathing crack. They observed that several sub or super- harmonic components are present in the spectra. However, the phenomena causing these several harmonics to emerge in the experiment have not been explained fully. The main reason for this is that the nonlinear vibration of a beam with a breathing crack is not fully understood, and that the inherent relation between the sub- and super-harmonic resonance responses and the parameters of the cracked beam require further investigation.

To fully understand the vibration characteristics caused by the breathing crack, increasing research activities have been focused on the use of nonlinear models for cracked beams. The modelling of the breathing crack is the key step that provides a solid foundation for the following analysis. There are two common closing crack models, namely a “square-type” model and a “simple periodic function” model [26]. The first model assumes that the breathing crack can only be fully open or fully closed, and the states of partial closure or opening are excluded [27,28,29]. For example, Chu and Shen [27] presented the square-type model to simulate the stiffness change, and they formulated a closed-form solution of the cracked beam. Bayly [30] used the perturbation approach to predict the amplitude of spectral peaks of the weakly bilinear oscillator. Chondrous et al. [31] developed a dynamic model of a cracked beam with a breathing crack by using consistent one-dimensional theory. Their results showed that the natural frequencies obtained from analysis are in accordance with the experimental results. However, this assumption is not accurate in most cases. For fatigue due to crack breathing, a partial crack closure occurs in reality. In particular, Clark et al. [32] and Rezaee et al. [33] have experimentally measured the state of partial closure. The second model simulates the stiffness changes as a simple trigonometric function [26] and has been adopted in [15,34]. The crack in this simplified model can be not only in a fully open state or fully closed state, but also in a state of being partially open or closed. However, this model can be applied only to the case of the fundamental mode [35].

The most common models for the vibration analysis of a cracked beam include the finite element model (FEM), in which a special cracked beam finite element is used. On that basis, several approaches in conjunction with the FEM have been proposed. The commonly-used approaches include the direct numerical integration method [36], the frequency response function (HRF) [20], the harmonic balance method [37], and mode superposition [38]. For example, Ruotolo et al. [20] used the FEM combined with the higher-order HRF to analyze the vibrational response of a cracked cantilever beam with a closing crack, and they found that the appearance of super-harmonic resonances strongly depends upon the crack size and position. Qian et al. [36] established the FEM of a cracked beam, and they compared the time response of the beam between the open crack and the breathing crack by numerical solution. Although the FEM of the cracked beam is by far the most commonly used method, this model has certain drawbacks which limit its application. For example, the FEM of a cracked beam with a breathing crack does not consider the states of partial crack closure. Besides this, the solving of the FEM of a beam with a breathing crack is relatively complex and time-consuming [25].

In this paper, a new stiffness model that considers the influence of partial crack closure is proposed to closely represent the stiffness variation of the cracked beam. Based on the proposed stiffness model, the equation of motion of a beam with a breathing crack subjected to a harmonic excitation is established. The multiple-scales method is applied to determine the resonance conditions. The relationships between the nonlinear behaviors and the system parameters are examined via both the numerical simulation and experimental study. The study shows that the frequent factors [39,40,41] may be used as a qualitative damage indicator for crack detection.

2. The Breathing Crack Model

A data set of load P and crack opening displacement (COD) u of a beam with a fatigue crack was experimentally measured using four-point bend testing by Clark et al. [32]. The COD gauge was fixed on the fatigue crack starter mouth of the specimen, which was used to measure the relative amount of the COD u, as shown in Figure 1. The strain gauge was attached to the uncracked side of the beam specimen, which was used to monitor the applied load. Based on the mechanics of materials, the relation between the amount of the COD u at y₁ and the strain at y₂ can be written as

$${\epsilon }_y=-\frac{y_{2}u}{y_1{u}_0}=-\frac{u}{u_0}.$$

(1)

where y₁ is the distance from the location of COD gauge to the neutral axis, y₂ is the distance from the location of the strain gauge to the neutral axis, and u₀ is the span of the COD gauge. Equation (1) indicates that the strain ε_y is proportional to the COD u. In other words, the curve P/u and P/ε_y will be similar. For convenience, a set of load P and strain ε_y data points is used to determine the breathing crack function.

As shown in Figure 2, the experimental data points between the load P and COD u are expressed by a “scatter” plot. These data are used to fit a three-section curve. The low section and upper section are a first-order polynomial (straight line), while the middle section is a second-order polynomial (quadratic) indicated by the red line in Figure 2. It is clear that the state of crack closing can be divided into three regions: fully closed, partially closed and fully open. This observation was also made in Cheng et al. [15] and Rezaee and Hassannejad [33]. From the fitting curves shown in Figure 2, it can be seen that the fully open and fully closed regions (states) are a straight line; i.e., the characteristics of the load P versus the COD u obey Hook’s law in these two regions. However, the state of the partially closed region does not exhibit such a linear characteristic, and a transition occurs if the crack begins to close up. When the COD is transferred from tensile state to compressive state, the crack surfaces experience a preliminary contact. As the COD reaches a certain compression value, the crack becomes fully closed. The boundary between the partially closed state and the fully closed state is named the “knot”, and the corresponding load P and COD u are denoted by P_k and u_k, respectively.

The stiffness of the cracked beam can be determined by finding the slope of the P–u curve at any displacement u [15,33]. Thus, the bending stiffness of a beam with a breathing crack k_t can be described in the form:

$$k_t=\{\begin{array}{l}{k}_{c}u<u_k,\\ c_{0}u+k_o{u}_k\le u<0,\\ k_{o}u\ge 0.\end{array}$$

(2)

where k_c and k_o are the bending stiffness of the cracked beam with a fully closed crack and fully open crack, respectively; c₀ is the unknown coefficient.

If the crack repetitively opens and closes during the vibration, the bending stiffness k_t will vary periodically, as shown in Figure 3. Therefore, the bending stiffness of a beam with a breathing crack can be written as the complex Fourier series in the form

$$k_t=k_o+f\left(t\right)(k_c-k_o).$$

(3)

where t is the time, f(t) is the breathing crack function which describes the process of crack opening and closing, and

$$f(t)=\frac{a_0}{2}+{\displaystyle \sum _{n=1}^N(a_n\cos n\omega t+b_n\cos n\omega t)}.$$

(4)

where ω is the frequency of the harmonic force acting on the cracked beam, a₀ is the undetermined constant, N is the number of harmonic terms used to model f(t), and a_n and b_n are the appropriate constants of the n-th order term, respectively. The constants a₀, a_n and $b_n$ can be determined by

$$\{\begin{array}{l}{a}_0=\frac{2}{T}{\displaystyle {\int }_0^{T}k(t)dt,}\\ a_n=\frac{2}{T}{\displaystyle {\int }_0^{T}k(t)\cos (n\omega )dt,}\\ b_n=\frac{2}{T}{\displaystyle {\int }_0^{T}k(t)\sin (n\omega )dt.}\end{array}$$

Using the trigonometric identity, Equation (4) can then be transformed into

$$f(t)=\frac{a_0}{2}+{\displaystyle \sum _{n=1}^N{k}_n\cos (n\omega t+{\phi }_n)}.$$

(5)

where

$$\{\begin{array}{l}{k}_n=\sqrt{a_n^2+b_n^2},\\ {\phi }_n=\arctan \frac{a_n}{b_n}.\end{array}$$

3. Stiffness Matrix of a Beam Element with a Breathing Crack

Figure 4 shows a beam element with a breathing crack. The element is defined by two nodes. Each node has three degrees of freedom (DOFs), where u₁ and u₄ are the longitudinal displacements of the node $i$ and the node $i+1$, respectively; u₂ and u₅ are the lateral displacement of the node $i$ and the node $i+1$, respectively; and u₃ and u₆ are the rotation angle of the node $i$ and the node $i+1$, respectively. The generalized coordinate of a beam element with a breathing crack u is expressed as

$$\mathbf{u}={\left\{u_1{u}_2{u}_3{u}_4{u}_5{u}_6\right\}}^{\mathrm{T}}.$$

According to the finite element theory, the transverse deflection $W(x)$ and longitudinal deflection $V(x)$ at an arbitrary neutral axis point of the cracked bam element can be expressed as [42]

$$\{\begin{array}{l}W(x)={\displaystyle \sum _i{B}_i(x)u_i,}i=2,3,5,6.\\ V(x)={\displaystyle \sum _j{B}_j(x)u_j,}j=1,4.\end{array}$$

(6)

where B₁(x), B₂(x), B₃(x), B₄(x), B₅(x), B₆(x) are the Hamiltonian interpolation functions, respectively, and

$$\{\begin{array}{l}{B}_1=1-\frac{x}{l},\\ B_2=1-\frac{3}{l^2}{x}^2+\frac{2}{l^3}{x}^3,\\ B_3=x-\frac{2}{l}{x}^2+\frac{1}{l^2}{x}^3,\\ B_4=\frac{x}{l},\\ B_5=\frac{3}{l^2}{x}^2-\frac{2}{l^3}{x}^3,\\ B_6=-\frac{1}{l}{x}^2+\frac{2}{l^2}{x}^3.\end{array}$$

where l is the length of the beam element with a breathing crack.

According to the mechanics of materials, the strain at an arbitrary neutral axis ${\epsilon }_x$ is composed of two parts, namely the bending strain ${\epsilon }_{{}_x}^b$ and the longitudinal strain ${\epsilon }_x^l$. Thus, the strain at $r$ shown in Figure 4 can be written as [42]

$${\epsilon }_r={\epsilon }_{{}_r}^l+{\epsilon }_{{}_r}^b=\frac{\partial V(x)}{\partial x}-r\frac{{\partial }^{2}W(x)}{\partial x^2}.$$

(7)

Applying Equation (6) to Equation (7) yields

$${\epsilon }_r=-\frac{1}{l}{u}_i+\frac{1}{l}{u}_{i+1}-r[(\frac{12x}{l^3}-\frac{6}{l^2})v_i+(\frac{6x}{l^2}-\frac{4}{l})v_{i+1}+(-\frac{12x}{l^3}+\frac{6}{l^2}){\theta }_i+(\frac{6x}{l^2}-\frac{2}{l}){\theta }_{i+1}].$$

(8)

For the sake of calculation and experimental measurement, the strain at r is used here to distinguish the state of the breathing crack. When the strain at r is in compression (ε_r ≤ 0), the crack is fully open. When the strain is transferred from compressive to tensile, the crack surfaces begin a preliminary contact. The crack surface acts as the partial crack closure state until the strain reaches the boundary value ε_k. The boundary strain value can be found from the intersection of the curve of the partially closed crack with the curve of the fully closed crack. When the strain ε_r is tensile and greater than ε_k (ε_r ≥ ε_k), the crack fully closes. Therefore, the stiffness matrix of the beam element with a breathing crack k_b can be expressed as

$${\mathbf{k}}_b=\{\begin{array}{l}{\mathbf{k}}_o{\epsilon }_r\le 0,\\ {\mathbf{k}}_o+\frac{{\epsilon }_r}{{\epsilon }_k}({\mathbf{k}}_c-{\mathbf{k}}_o){\epsilon }_k>{\epsilon }_r>0,\\ {\mathbf{k}}_c{\epsilon }_k\le {\epsilon }_r.\end{array}$$

(9)

where k_c is the stiffness matrix of the cracked beam element when the crack is fully closed, or the stiffness matrix of the beam element without a crack; k_o is the stiffness matrix of the cracked beam element when the crack fully opens, which can be expressed as [26]

k_o=Φ (c₀+ c₁)⁻¹Φ ^T.

(10)

where c₀ is the flexibility matrix of the beam element without a crack, c₁ is the additional flexibility matrix caused by the crack, the subscript t is the transformation of a matrix, and Φ is the transformation matrix given as

$$\mathsf{\Phi }={\left[\begin{array}{cccccc}0& 0& -1& 0& 0& 1\\ 0& -1& -l& 0& 1& 0\\ -1& 0& 0& 1& 0& 0\end{array}\right]}^{\mathrm{T}}$$

When the cracked beam vibrates repetitively during vibration, the stiffness matrix k_b can be written in the Fourier series:

k_b = k_o + f(t)(k_c−k_o).

(11)

4. Finite Element Model of a Beam with a Breathing Crack

The cracked beam is discretized by one cracked beam element and (n−1) non-crack beam elements, as shown in Figure 5, where n is the number of elements. Each of the non-crack elements has two nodes and three DOFs at each of the nodes, and the cracked beam element has the same node numbers and DOF numbers as the non-crack elements. Based on the assumption that the changes of damping caused by the crack are relatively small compared to the total damping of the cracked beam, the changes of the damping can be neglected. In addition, it is also assumed that the crack has no effect on the mass of the system. Thus, the nonlinear equation of motion of the cracked beam with a breathing crack is expressed as

$$\mathbf{M}\ddot{\mathbf{U}}+\mathbf{C}\dot{\mathbf{U}}+[{\mathbf{K}}_{01}-{\mathbf{K}}_{02}f(t)]\mathbf{U}=\mathbf{F}.$$

(12)

where M is the mass matrix of the cracked beam, C is the damping matrix of the cracked beam, K₀₁ is the stiffness matrix of the undamaged beam, and K₀₂ is the stiffness matrix caused by the crack; U is the generalized coordinate vector of the cracked beam, and F is the applied force vector.

Substituting Equation (5) into Equation (12), the dynamic equation of the cracked beam can be written as

$$\mathbf{M}\ddot{\mathbf{U}}+\mathbf{C}\dot{\mathbf{U}}+({\mathbf{K}}_{001}-{\mathbf{K}}_t)\mathbf{U}=\mathbf{F}.$$

(13)

where

$${\mathbf{K}}_{001}={\mathbf{K}}_{01}-\frac{a_0{\mathbf{K}}_{02}}{2},{\mathbf{K}}_t={\mathbf{K}}_{02}{\displaystyle \sum _{n=1}^N{k}_n\cos (n\omega t+{\phi }_n)}.$$

Equation (13) shows that the stiffness matrix K_t is time-variant due to the repetitive crack opening and closing. Thus, Equation (13) can be regarded as a linear time-varying equation.

5. Analysis of a Beam with a Breathing Crack

The resonance conditions of the cracked beam are determined by using the multiple-scale method. This method assumes that the system has a weak damping and a weak exciting force. To facilitate the analysis, a parameter $\epsilon$ is introduced as the designator for the damping matrix C, stiffness matrix K_t and the force vector F. Now, Equation (13) can be rewritten as

$$\mathbf{M}\ddot{\mathbf{U}}+{\mathbf{K}}_{001}\mathbf{U}=\epsilon \overline{\mathbf{F}}-\epsilon \overline{\mathbf{C}}\dot{\mathbf{U}}+\epsilon {\mathbf{K}}_{02}\mathbf{U}{\displaystyle \sum _{n=1}^N{\overline{k}}_n\cos (n\omega t+{\phi }_n}).$$

(14)

where

$$\overline{\mathbf{F}}=\frac{\mathbf{F}}{\epsilon },\overline{\mathbf{C}}=\frac{\mathbf{C}}{\epsilon },{\overline{k}}_n=\frac{k_n}{\epsilon }.$$

To predict the steady state response of the cracked beam, a mode superposition method can be used. The linear transformation is introduced as

$$\{\begin{array}{l}\mathbf{U}=\mathsf{\psi }\mathsf{\eta },\\ \dot{\mathbf{U}}=\mathsf{\psi }\dot{\mathsf{\eta }},\\ \ddot{\mathbf{U}}=\mathsf{\psi }\ddot{\mathsf{\eta }}.\end{array}$$

(15)

where ψ is the matrix of the modal transformation that can be obtained from the FEM of the cracked beam, and η is the modal coordinate vector.

Substituting $\mathbf{U}$, $\dot{\mathbf{U}}$, and $\ddot{\mathbf{U}}$ given in Equation (15) into Equation (14) and premultiplying it by ${\mathsf{\psi }}^{\mathrm{T}}$ yields

$$\ddot{\mathsf{\eta }}+{\mathbf{K}}_{0001}\mathsf{\eta }=\epsilon {\mathbf{F}}_0-\epsilon {\mathbf{C}}_0\dot{\mathsf{\eta }}+\epsilon {\mathbf{K}}_{002}\mathsf{\eta }{\displaystyle \sum _{n=1}^N{\overline{k}}_n\cos (n\omega t+{\phi }_n)}.$$

(16)

where

$${\mathbf{K}}_{0001}={\mathsf{\psi }}^{\mathrm{T}}{\mathbf{K}}_{001}\mathsf{\psi },{\mathbf{C}}_0={\mathsf{\psi }}^{\mathrm{T}}\overline{\mathbf{C}}\psi ,{\mathbf{F}}_0={\mathsf{\psi }}^{\mathrm{T}}\overline{\mathbf{F}},{\mathbf{K}}_{002}={\mathsf{\psi }}^{\mathrm{T}}{\mathbf{K}}_{02}\mathsf{\psi }.$$

The solution of Equation (16) can be expressed as

$$\mathsf{\eta }={\mathsf{\eta }}_0(T_0,T_1,\dots ,T_m)+\epsilon {\mathsf{\eta }}_1(T_0,T_1,\dots ,T_m)+\cdots +{\epsilon }^m{\mathsf{\eta }}_m(T_0,T_1,\dots ,T_m).$$

(17)

where the independent time variable is defined as [43]

$$T_m={\epsilon }^m\tau ,m=0,1,2,\dots ,$$

(18)

where $T_0$ is the fast timescale, and $T_1$ and $T_2$ are slow timescales that are much slower than $T_0$; τ is the normalized time.

According to the multiple-scale method, the time derivatives with respect to τ are given by [43]

$$\{\begin{array}{l}\frac{d}{d\tau }=D_0+\epsilon D_1+{\epsilon }^2{D}_2+\cdots ,\\ \frac{d^2}{d{\tau }^2}=D_0^2+2\epsilon D_0{D}_1+{\epsilon }^2(D_1^2+2D_0{D}_2)+\cdots .\end{array}$$

(19)

where $D_m=\frac{\partial }{\partial T_m}$, is partial derivative operator symbol.

The first-order approximation of η can be written as

$$\mathsf{\eta }={\mathsf{\eta }}_0(T_0,T_1)+\epsilon {\mathsf{\eta }}_1(T_0,T_1).$$

(20)

The r-th element in the vector η can be expressed as

η_r = η_0r + εη_1r, r = 1, 2, …, q.

(21)

where q is the number of DOFs.

After substituting Equations (19) and (20) into Equation (16), the following hierarchy of equations can be obtained for ε⁰ and ε¹, respectively:

$$D_0^2{\mathsf{\eta }}_0+{\mathbf{K}}_{0001}{\mathsf{\eta }}_0=0$$

(22)

$$D_0^2{\ddot{\mathsf{\eta }}}_1+{\mathbf{K}}_{0001}{\mathsf{\eta }}_1=-2D_0{D}_1{\mathsf{\eta }}_0+{\mathbf{F}}_0-C_0{D}_0{\mathsf{\eta }}_0+{\displaystyle \sum _{n=1}^N{\mathbf{K}}_{002}{\mathsf{\eta }}_0{\overline{k}}_n\cos (n\omega t+{\phi }_n)}$$

(23)

Then, the solution to Equation (22) can be written as

$${\eta }_0{}_r=A_r(T_1)e^{j{\omega }_r{T}_0}+{\overline{A}}_r(T_1)e^{-j{\omega }_r{T}_0}$$

(24)

where $A_r(T_1)$ is a complex valued amplitude of the r-th order of the cracked beam, and ${\overline{A}}_r(T_1)$ indicates a complex conjugate of $A_r(T_1)$; ω_r is the r-th order natural frequency, which can be determined by the mass matrix M the stiffness matrix of the uncracked beam K₀₁, and the stiffness matrix K₀₂ of Equation (12).

Substituting Equation (24) into Equation (23) yields

$$\begin{array}{ll}{D}_0^2{\eta }_{r1}+{\omega }_r^2{\eta }_{r1}=& -2j{\omega }_r({A^{\prime }}_r+c_{0r}{A}_r)e^{j{\omega }_r{T}_0}+{\displaystyle \sum _{n=1}^N{\displaystyle \sum _{p=1}^q\frac{{\overline{k}}_n}{2}{({\mathbf{K}}_{002})}_{p,r}{A}_r{e}^{j{\omega }_r{T}_0}{e}^{j(n\omega T_0+{\phi }_n)}}}+\\ & {\displaystyle \sum _{n=1}^N{\displaystyle \sum _{p=1}^q\frac{{\overline{k}}_n}{2}{({\mathbf{K}}_{002})}_{p,r}{\overline{A}}_r{e}^{-j{\omega }_r{T}_0}{e}^{j(n\omega T_0+{\phi }_n)}}}+({\mathbf{F}}_0{)}_r{e}^{j\omega T_0}+cc.\end{array}$$

(25)

where the prime represents a time derivative of slow timescale T₁, and cc. represents the complex conjugate of the proceeding terms.

According to the definition by Wang et al. [39] and Li et al. [40,41], the frequency of each term in the first-order approximation solution of a nonlinear system is defined as the frequent factor. Therefore, the corresponding frequent factors of Equation (25) are as follows:

ω, nω ± ω_r (r = 1, 2, …, q; n = 1, 2, …, N),

It is seen that the steady state responses consist of several sub- and super-harmonic components, whose frequencies related to the natural frequencies are nω ± ω_r. Based on the theory of nonlinear resonances [44], the cracked beam may be at resonance under the following conditions:

(1)

If ω ≈ ω_r, primary resonance will occur.

(2)

If 6ω − ω_r ≈ ω_r or ω ≈ 1/3 ω_r, third-order super-harmonic resonance will occur.

(3)

If 4ω − ω_r ≈ ω_r or ω ≈ 1/2ω_r, second-order super-harmonic resonance will occur.

(4)

If ω − ω_r ≈ ω_r or ω ≈ 2 ω_r, 1/2th-order sub-harmonic resonance will occur.

(5)

If nω ± ω_r ≈ ω_i (r ≠ i), combination resonance will occur.

(6)

If two types of the resonance occur simultaneously, multiple resonance will occur.

As shown above, the natural frequency ω_r is determined by both the structural parameters and the crack parameters. Apparently, the frequent factors are related to structural parameters, the crack parameters and the exciting force. This indicates that, given an exciting force, the nonlinear vibration of the cracked beam is influenced not only by the parameters of the uncracked beam but also by the crack parameters.

6. Solution of the Nonlinear Equation

In general, the analysis accuracy can be satisfied by considering the first few modes of the system [42]. In order to improve the solving efficiency, the component mode synthesis (CMS) method [45] is used in this paper. Using the truncated normal mode matrix ψ_C to reduce the number of DOFs, the reduced order model of the beam with a breathing crack can be written as

$${\mathbf{M}}_C{\ddot{\mathsf{\eta }}}_C+{\mathbf{C}}_C{\dot{\mathsf{\eta }}}_C+{\mathbf{K}}_C{\mathsf{\eta }}_C-{\mathbf{K}}_T{\mathsf{\eta }}_C{\displaystyle \sum _{n=1}^N{k}_n\cos (n\omega t+{\phi }_n)}={\mathbf{F}}_C$$

(26)

where ${\mathbf{M}}_C={\mathsf{\psi }}_C^{\mathrm{T}}\mathbf{M}{\mathsf{\psi }}_C$, ${\mathbf{C}}_C={\mathsf{\psi }}_C^{\mathrm{T}}\mathbf{C}{\mathsf{\psi }}_C$, ${\mathbf{K}}_C={\mathsf{\psi }}_C^{\mathrm{T}}{\mathbf{K}}_{001}{\mathsf{\psi }}_C$, ${\mathbf{K}}_T={\mathsf{\psi }}_C^{\mathrm{T}}{\mathbf{K}}_{02}{\mathsf{\psi }}_C$, and ${\mathbf{F}}_C={\mathsf{\psi }}_C^{\mathrm{T}}\mathbf{F}$,$\mathbf{U}={\mathsf{\psi }}_C{\mathsf{\eta }}_C$.

It should be noted that the CMS method is usually employed to treat larger structures as an assembly of multiple substructures [45]. It is used here because it provides a convenient modeling framework to reduce the DOF number of the FEM of the cracked beam.

7. Results and Comparisons

To validate the analytical results, the forced responses of a beam with a breathing crack are studied by both experiment and numerical calculation. Besides this, to evaluate the accuracy of the proposed model, the experimental results and numerical results are compared.

7.1. Experimental Platform

Fatigue crack preparation is a key procedure for the following vibration testing. First, a beam with a notch is cut from a C45E steel plate, and its dimensions in mm are given in Figure 6a. The elasticity modulus of the material is $E$ = 206 GPa and the mass density is $\rho$ = 7850 kg/m³. The actual fatigue crack is created by a servo hydraulic machine (SDS 200), as shown in Figure 6b. The bending fatigue testing for crack growth is conducted in accordance with the specification ASTM-E647. The constant amplitude loading test is performed with a loading frequency of 20 Hz and a stress ratio of 0.1. The crack length is 3.825 mm, measured by an optical microscope (AO-AF216C), as shown in Figure 6c. Then, to remove the notch starter, the height of the specimen is machined to 12 mm by the electrical discharge machine (EDM).

Figure 7 shows the experimental setup. It consists of three parts: the exciting shaker system, the vibrating system, and the data acquisition system. The vibrating system consists of the cracked beam specimen, the solid base, the force sensor (VL25), the threaded stem, the shaker (HEA-50) and the aluminum block, as shown in Figure 7a. The cracked beam specimen is clamped on the solid base by bolts. A force sensor weighing 0.05 kg is fixed at the free end of the beam by bolts. An accelerometer weighing 0.01 kg is attached at the free end of the cantilever beam. The harmonic force generated by the shaker is transmitted to the beam through a threaded stem which connects the force sensor and the shaker by a screw. The shaker system consists of the signal generator, the power amplifier (HEAS-50) and the shaker. The data acquisition system consists of an accelerometer (JF-2050) and a PC equipped with a data acquisition board (DEWE-201), as shown in Figure 7b.

7.2. Experimental Results

The dimensions in mm of the cracked cantilever beam are shown in Figure 8. The distance from the crack location to the fixed end is 176 mm. The crack depth and width are 3.8 mm and 12 mm, respectively. From the free vibrations of the cracked beam excited by a hammer impact, the first natural frequency f₁ is found to be 57.4 Hz.

The forced vibration consists of three groups. The external harmonic force of each test is 20 N, and the excited frequency ratios λ = f_ω/f₁ are considered to be λ = 1/3, 1/2, 2. The sampling frequency is set as 10,000 Hz. The spectrum domain of the acceleration signals for the three cases λ = 1/3, 1/2, and 2 are depicted in Figure 9, Figure 10 and Figure 11, respectively. The spectrum analysis is conducted by fast Fourier transform (FFT). To evaluate the nonlinear characteristics, a relative harmonic amplitude is introduced [22]:

$$h_r=\frac{H_s}{H_f}$$

(27)

where H_s is the magnitude of the sub-harmonic or super-harmonic component, and H_f is the FFT magnitude of the harmonic component corresponding to the exciting frequency.

Figure 9 shows that when f_ω = 1/3f₁ = 19.1 Hz, the spectrum exhibits rich harmonic contents. The harmonic at f₁ = 3f_ω = has the largest amplitude, and the corresponding value of h_r is 132%. The magnitude of the harmonic at the exciting frequency f_ω (1/3f₁) is slightly smaller than that of the third harmonic component. Obviously, the several components at the frequency 2f_ω (f ₁ − f_ω), 3f_ω (f ₁), 4f_ω (f₁ + f_ω), and 5f_ω (f₁ + 2f_ω) coincide with the frequent factors nω ± ω₁ (n = 1, 2), which validate the correctness of the analytical results from the multiple-scale method.

Figure 10 shows that when f_ω = 1/2f₁ = 28.7 Hz, the response also contain several harmonics. The amplitude reaches its maximum at the second harmonic 2f_ω (f₁), which is higher than the harmonic component at the exciting frequency f_ω (1/2f₁). In this case, the corresponding h_r is 179%. Besides this, the harmonics 3f_ω (f₁ + f_ω) and 4f_ω (f₁ + 2f_ω) are also observed. It should be noticed that these harmonic components coincide with the frequent factors nω ± ω₁ (n = 1, 2).

Figure 11 shows that when f_ω = 2f₁ = 114.8 Hz, the response has only two dominant harmonics in the frequency range of consideration. The peak at the sub-harmonic with f₁ (1/2f_ω) is smaller than the excitation harmonic component at the exciting frequency f_ω. The corresponding h_r is 23%. The experimental results show that the sub-harmonic resonance occurs at f_ω = 2f₁.

7.3. Numerical Analysis

To validate the proposed model, the forced response obtained from the proposed model are compared to the experiment results. The structure of the experimental beam sketched in Figure 8 are also discretized in seven beam elements with three DOFs. The lengths of the beam elements are, l₁ = 40 mm, l₂ = 40 mm, l₃ = 60 mm, l₄ = 70 mm, l₅ = 70 mm, l₆ = 70 mm, and l₇ = 56 mm, respectively. The breathing crack is set on element 4. The damping is assumed to be a proportional model with a damping ratio of ξ = 0.05, which is estimated by applying the half-power bandwidth method to the forced response. The boundary strain is ε_k = 0.004. The forced response is obtained by solving Equation (26) through the MATLAB program by adopting the fourth-order Runge–Kutta method. The forced responses which are obtained from the proposed model for different excited frequencies are plotted with a dashed line in Figure 9, Figure 10 and Figure 11. It can be seen that the numerical results agree well with the experimental results. The numerical results have the same harmonics as those in the experimental results. The harmonic component frequencies obtained from the numerical results also coincide with the analyzed frequent factors, which validates the correctness of the analytical results.

To determine the differences between the force responses from the proposed model and those from the experiment,

Table 1. Comparison of the forced response of dominant harmonic amplitude from the proposed method and experiment for different driving frequencies.

Crack		Driving Frequency	Proposed Method	Experiment	Difference (%)
X/L	a/h	Ratio λ
		1/3	130	132	1.52
0.37	0.32	1/2	188	179	5.03
		2	22	23	4.35

compares the relative amplitude of the dominant harmonic h_r obtained from proposed model and those from the experimental results. For the three forcing frequencies f_w = 1/3f₁, 1/2 f₁, and 2 f₁, the relative amplitude of the dominant harmonic components is largest when the forcing frequency is f_w = 1/2 f₁. This indicates that the most effective forcing frequency for crack detection is f_ω = 1/2f₁. As shown in Table 1, the differences between the relative harmonic amplitude determined by the proposed method and those from the experiment range from 1.52% to 5.03%, and the maximum relative error is no larger than 6%. Such a good agreement validates the accuracy of the proposed model.

8. Conclusions

In this paper, we have developed a new stiffness model that considers the influence of partial crack closure. Using this model, we have derived a nonlinear dynamic equation of a beam with a breathing crack by the finite element method. Then, we applied the multiple-scales method to analyze the approximate steady state response of the cracked beam to a harmonic excitation. We have shown that the steady state responses consist of multiple harmonics. Referred to as frequent factors, the frequencies of these harmonic terms can be used to reveal whether a beam has a crack or not. The main findings of the study can be summarized as follows:

• The study shows that the nonlinear vibration of a cracked beam is influenced not only by the parameters of its structure, but also by the crack parameters and excited force parameters. Under some circumstances, this system presents third-order super-harmonic resonance, second-order super-harmonic resonance, 1/2th-order sub-harmonic resonance, combination resonance, and multiple resonance.
• The frequent factors show that several sub- and super-harmonic components exist in the dynamic response, whose frequencies related to the natural frequencies are nω ± ω_r (r = 1, 2, …, q; n = 1, 2, …, N). The analytical results agree well with the experimental data and numerical results. Therefore, the frequent factor can be used as a qualitative damage indicator for crack detection.
• When the beam is excited at one of those frequency factors, the spectrum exhibits obvious nonlinear characteristics. In particular, when the cracked beam is excited by the forcing frequency f_ω = 1/2f₁, the nonlinear behavior is most visible, indicating that this frequency is the most efficient one to reveal damage.
• The relative harmonic amplitudes calculated by the proposed method agree well with the experimental values for the three cases under consideration. This shows that the proposed stiffness model is reliable.