Nonlinear Vibration Analysis of a Beam with a Breathing Crack

: 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 sti ﬀ ness model that considers the inﬂuence of the partial crack closure is proposed to model the sti ﬀ ness variation of the cracked beam. Based on the ﬁnite 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 a ﬀ ected by both the structural parameters and the crack parameters.


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.

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 1 and the strain at y 2 can be written as where y 1 is the distance from the location of COD gauge to the neutral axis, y 2 is the distance from the location of the strain gauge to the neutral axis, and u 0 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.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 3 of 17 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 timeconsuming [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.

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 y1 and the strain at y2 can be written as where y1 is the distance from the location of COD gauge to the neutral axis, y2 is the distance from the location of the strain gauge to the neutral axis, and u0 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 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. 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 Pk and uk, respectively.

Partially closed
Fully closed

Fully open
Crack open displacement u (mm) Load P (kN) (uk, pk) knot 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 kt can be described in the form: where kc and ko are the bending stiffness of the cracked beam with a fully closed crack and fully open crack, respectively; c0 is the unknown coefficient. If the crack repetitively opens and closes during the vibration, the bending stiffness kt 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 where t is the time, f(t) is the breathing crack function which describes the process of crack opening and closing, and 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: 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 0 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 where t is the time, f (t) is the breathing crack function which describes the process of crack opening and closing, and (a n cos nωt + b n cos nωt).
where ω is the frequency of the harmonic force acting on the cracked beam, a 0 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 0 , a n and b n can be determined by Using the trigonometric identity, Equation (4) can then be transformed into where where ω is the frequency of the harmonic force acting on the cracked beam, a0 is the undetermined constant, N is the number of harmonic terms used to model f(t), and an and bn are the appropriate constants of the n-th order term, respectively. The constants a0, an and can be determined by Using the trigonometric identity, Equation (4) can then be transformed into  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 u1 and u4 are the longitudinal displacements of the node i and the node i + 1 , respectively; u2 and u5 are the lateral displacement of the node i and the node i + 1 , respectively; and u3 and u6 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   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 1 and u 4 are the longitudinal displacements of the node i and the node i + 1, respectively; u 2 and u 5 are the lateral displacement of the node i and the node i + 1, respectively; and u 3 and u 6 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 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]

Stiffness Matrix of a Beam Element with a Breathing Crack
where are the Hamiltonian interpolation functions, respectively, and where l is the length of the beam element with a breathing crack.
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] where B6(x) are the Hamiltonian interpolation functions, respectively, where l is the length of the beam element with a breathing crack. ε . Thus, the strain at r shown in Figure 4 can be written as [42] ( ) ( ) Applying Equation (6) to Equation (7) yields 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 According to the mechanics of materials, the strain at an arbitrary neutral axis ε x is composed of two parts, namely the bending strain ε b x and the longitudinal strain ε l x . Thus, the strain at r shown in Figure 4 can be written as [42] Applying Equation (6) to Equation (7) yields 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 Appl. Sci. 2019, 9, 3874 7 of 17 (ε r ≥ ε k ), the crack fully closes. Therefore, the stiffness matrix of the beam element with a breathing crack k b can be expressed as 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] where c 0 is the flexibility matrix of the beam element without a crack, c 1 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 When the cracked beam vibrates repetitively during vibration, the stiffness matrix k b can be written in the Fourier series:

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 M . .
where M is the mass matrix of the cracked beam, C is the damping matrix of the cracked beam, K 01 is the stiffness matrix of the undamaged beam, and K 02 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.
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 kb can be expressed as where kc 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; ko is the stiffness matrix of the cracked beam element when the crack fully opens, which can be expressed as [26] ko=Φ (c0+ c1) -1 Φ T .
where c0 is the flexibility matrix of the beam element without a crack, c1 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 0 0 1 0 0 1 When the cracked beam vibrates repetitively during vibration, the stiffness matrix kb can be written in the Fourier series: kb = ko +f(t)(kc-ko). (11)

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 where M is the mass matrix of the cracked beam, C is the damping matrix of the cracked beam, K01 is the stiffness matrix of the undamaged beam, and K02 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. where k n cos(nωt + ϕ 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.

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 ε 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 where To predict the steady state response of the cracked beam, a mode superposition method can be used. The linear transformation is introduced as 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.

U given in Equation
The solution of Equation (16) can be expressed as where the independent time variable is defined as [43] T m = ε m τ, m = 0, 1, 2, . . . , 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] where , is partial derivative operator symbol.
Appl. Sci. 2019, 9, 3874 9 of 17 The first-order approximation of η can be written as The r-th element in the vector η can be expressed as 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 ε 0 and ε 1 , respectively: Then, the solution to Equation (22) can be written as where A r (T 1 ) is a complex valued amplitude of the r-th order of the cracked beam, and 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 01 , and the stiffness matrix K 02 of Equation (12). Substituting Equation (24) into Equation (23) yields k n 2 (K 002 ) p,r A r e jω r T 0 e j(nωT 0 +ϕ n ) + N n=1 q p=1 k n 2 (K 002 ) p,r A r e −jω r T 0 e j(nωT 0 +ϕ n ) + (F 0 ) r e jωT 0 + cc.
where the prime represents a time derivative of slow timescale T 1 , 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.
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.

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 where 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.

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.

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 ρ = 7850 kg/m 3 . 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.

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 1 is found to be 57.4 Hz.

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 f1 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ω/f1 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 Figures 9-11, respectively. The spectrum analysis is conducted by fast Fourier transform (FFT). To evaluate the nonlinear characteristics, a relative harmonic amplitude is introduced [22]: where Hs is the magnitude of the sub-harmonic or super-harmonic component, and Hf is the FFT magnitude of the harmonic component corresponding to the exciting frequency.  The forced vibration consists of three groups. The external harmonic force of each test is 20 N, and the excited frequency ratios λ = f ω /f 1 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 Figures 9-11, respectively. The spectrum analysis is conducted by fast Fourier transform (FFT). To evaluate the nonlinear characteristics, a relative harmonic amplitude is introduced [22]: 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.

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 f1 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ω/f1 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 Figures 9-11, respectively. The spectrum analysis is conducted by fast Fourier transform (FFT). To evaluate the nonlinear characteristics, a relative harmonic amplitude is introduced [22]: where Hs is the magnitude of the sub-harmonic or super-harmonic component, and Hf is the FFT magnitude of the harmonic component corresponding to the exciting frequency.    Figure 9 shows that when fω = 1/3f1 = 19.1 Hz, the spectrum exhibits rich harmonic contents. The harmonic at f1 = 3fω = has the largest amplitude, and the corresponding value of hr is 132%. The magnitude of the harmonic at the exciting frequency fω (1/3f1) is slightly smaller than that of the third harmonic component. Obviously, the several components at the frequency 2fω (f 1 − fω), 3fω (f 1), 4fω (f1 + fω), and 5fω (f1 + 2fω) coincide with the frequent factors nω ± ω1 (n = 1, 2), which validate the correctness of the analytical results from the multiple-scale method. Figure 10 shows that when fω = 1/2f1 = 28.7 Hz, the response also contain several harmonics. The amplitude reaches its maximum at the second harmonic 2fω (f1), which is higher than the harmonic component at the exciting frequency fω (1/2f1). In this case, the corresponding hr is 179%. Besides this,   Figure 9 shows that when fω = 1/3f1 = 19.1 Hz, the spectrum exhibits rich harmonic contents. The harmonic at f1 = 3fω = has the largest amplitude, and the corresponding value of hr is 132%. The magnitude of the harmonic at the exciting frequency fω (1/3f1) is slightly smaller than that of the third harmonic component. Obviously, the several components at the frequency 2fω (f 1 − fω), 3fω (f 1), 4fω (f1 + fω), and 5fω (f1 + 2fω) coincide with the frequent factors nω ± ω1 (n = 1, 2), which validate the correctness of the analytical results from the multiple-scale method. Figure 10 shows that when fω = 1/2f1 = 28.7 Hz, the response also contain several harmonics. The amplitude reaches its maximum at the second harmonic 2fω (f1), which is higher than the harmonic component at the exciting frequency fω (1/2f1). In this case, the corresponding hr is 179%. Besides this, Figure 11. FFT spectrum of the transverse acceleration at the tip of the cracked beam with λ = 2. Figure 9 shows that when f ω = 1/3f 1 = 19.1 Hz, the spectrum exhibits rich harmonic contents.
The harmonic at f 1 = 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 1 ) is slightly smaller than that of the third harmonic component. Obviously, the several components at the frequency 2f ω (f 1 − f ω ), 3f ω (f 1 ), 4f ω (f 1 + f ω ), and 5f ω (f 1 + 2f ω ) coincide with the frequent factors nω ± ω 1 (n = 1, 2), which validate the correctness of the analytical results from the multiple-scale method. Figure 10 shows that when f ω = 1/2f 1 = 28.7 Hz, the response also contain several harmonics. The amplitude reaches its maximum at the second harmonic 2f ω (f 1 ), which is higher than the harmonic component at the exciting frequency f ω (1/2f 1 ). In this case, the corresponding h r is 179%. Besides this, the harmonics 3f ω (f 1 + f ω ) and 4f ω (f 1 + 2f ω ) are also observed. It should be noticed that these harmonic components coincide with the frequent factors nω ± ω 1 (n = 1, 2). Figure 11 shows that when f ω = 2f 1 = 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 (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 1 .

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 1 = 40 mm, l 2 = 40 mm, l 3 = 60 mm, l 4 = 70 mm, l 5 = 70 mm, l 6 = 70 mm, and l 7 = 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 Figures 9-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 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 , 1/2 f 1 , and 2 f 1 , the relative amplitude of the dominant harmonic components is largest when the forcing frequency is f w = 1/2 f 1 . This indicates that the most effective forcing frequency for crack detection is f ω = 1/2f 1 . 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.

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 1 , 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.