Use of Bispectrum Analysis to Inspect the Non-Linear Dynamic Characteristics of Beam-Type Structures Containing a Breathing Crack

A breathing crack is a typical form of structural damage attributed to long-term dynamic loads acting on engineering structures. Traditional linear damage identification methods suffer from the loss of valuable information when structural responses are essentially non-linear. To deal with this issue, bispectrum analysis is employed to study the non-linear dynamic characteristics of a beam structure containing a breathing crack, from the perspective of numerical simulation and experimental validation. A finite element model of a cantilever beam is built with contact elements to simulate a breathing crack. The effects of crack depth and location, excitation frequency and magnitude, and measurement noise on the non-linear behavior of the beam are studied systematically. The result demonstrates that bispectral analysis can effectively identify non-linear damage in different states with strong noise immunity. Compared with existing methods, the bispectral non-linear analysis can efficiently extract non-linear features of a breathing crack, and it can overcome the limitations of existing linear damage detection methods used for non-linear damage detection. This study’s outcome provides a theoretical basis and a paradigm for damage identification in cracked structures.


Introduction
Preventing damage from occurring in engineering structures is difficult, due to the long-term effects of loads and environmental factors [1]. A typical example of damage is a crack. This not only affects the integrity of an engineering structure, but also poses a serious threat to its safety and service life [2]. In particular, the presence of cracks in beam structures is found in numerous cases, making basic research into the mechanical behavior of such beams [3] of great significance for diagnosing any structural damage and developing health-monitoring technology.
Traditional damage-detection methods based on signal-processing techniques, typically Fourier transform and power-spectrum analysis, are mostly applied by assuming that engineering structures behave in a linear manner [4]. However, a vast amount of damage that occur in actual engineering structures cause non-linear dynamic behaviors [5][6][7][8]. The non-linear dynamic behaviors are mainly represented by two types: (i) progressive degradation or deterioration of structural properties with time evolution. The yielding of frame structures under seismic action is an example of this type of non-linear damage [9][10][11]; (ii) generation of secondary vibrational resource complicating structural dynamic behavior. This type of non-linear damage is exemplified by breathing cracks [12][13][14] and material stratification [15,16]. Wavelet-based time-frequency methods are representative techniques for characterizing structural progressive degradation or deterioration, the first type of non-linear damage, with its ability to detect the relationship between frequencies and its time dependence, and has a certain noise robustness [17][18][19][20][21][22]. Nevertheless, the method suitable for characterizing the second type of non-linear damage, particularly the breathing crack, is pretty absent while desirable.
The essence of traditional power-spectrum analysis [23], which plays an important role in signal processing, is the analysis of second-order signal statistics. Power-spectrum analysis can only provide information about amplitude-frequency characteristics of signals, but cannot provide high-order statistical information, nor can it effectively suppress the influence of noise. Unlike power-spectrum analysis, bispectral analysis has been explored in applications of non-linear analysis. Collis et al. [24] proved that certain information, unavailable through the examination of second-order statistics such as the spectrum or correlation function, can be obtained by higher-order spectra.
Li et al. [25] applied bispectral analysis for the diagnosis of rotor crack. Yiakopoulos et al. [26] showed that bispectral analysis could effectively extract rolling bearings' groundfault characteristics.
In general, bispectrum-based damage diagnosis relies on the correlation between spectral peaks and damage types. Commonly used are the horizontal, diagonal, and central frequency sections of a slice spectrum [27][28][29], and the cepstrum, after the cluster lines are simplified to a single spectral line [30,31].
Moreover, some studies have proved that damage types can be identified according to the spectrum-peak frequency of the bispectrum, and that the degree of damage can be identified by texture features [32]. In existing research and engineering applications, the bispectrum is mostly used in damage detection in rotating machinery such as bearings [33,34], but has rarely been used for damage detection in civil engineering structures. Thus there is a strong incentive to explore the advantages of the bispectrum in non-linear damage detection in civil engineering structures.
For this paper, a finite element model of a beam with a breathing crack was built, using contact elements to simulate the opening and closing characteristics of the breathing crack. Subsequently, the non-linear characteristics of the breathing cracked beam under harmonic load, the influence of excitation frequency and magnitude, crack position and depth, and noise on the non-linear characteristics of the structure were studied systematically.
The results show that the bispectral analysis method demonstrates sensitivity to the non-linear features of the cracked beam, and can effectively overcome the limitations of existing linear damage-detection methods in non-linear damage detection.

Fundamentals
The effectiveness of traditional damage-detection methods, based on linear hypotheses, is largely limited in non-linear structures. With its capacity to tackle non-linear problems, bispectrum analysis has the lowest order among the high-order spectrum-analysis approaches. Its development is based on the power-spectrum-analysis method. Compared with linear damage-detection methods, bispectrum analysis contains richer damage information and is an effective tool for dealing with non-linear signals.

Definition of Bispectrum
It is assumed that the mean value of a real random signal is zero and its k-order cumulant c kx (τ 1 , τ 2 , . . . , τ k−1 ) is stationary. The k-order spectrum is defined as the (k-1)dimensional Fourier transform of c kx (τ 1 , τ 2 , . . . , τ k−1 ), expressed as The third-order spectrum, i.e., S 3x with k = 3, is also called the bispectrum, expressed as B x (ω 1 , ω 2 ). The estimated bispectrum [35,36] can be written as is the Fourier transform determined by the principal value sequence of {x(n)}; * represents the complex conjugate. B(ω 1 , ω 2 ) is a function with two independent frequencies as the independent variables, i.e., ω 1 and ω 2 , which contain the phase information. Since the bispectrum is symmetric, only the function parts within the regions of ω 1 ≤ ω 2 , 0 ≤ ω 2 ≤ f 2 /2 and ω 1 + 2ω 2 ≤ f s need to be calculated, where f s represents the sampling frequency.
In the following study, the bispectrum is treated by logarithm and normalization in order to amplify the harmonic components.

Identification of Non-Linear Systems
While the power spectrum represents the distribution of signal energy against frequency, the physical meaning of bispectrum is not explicit. Reference [37] explains the physical meaning of bispectrum: since the secondary moment of zero delay is the variance of the signal, and the third moment of zero delay is the skewness of the signal, the power spectrum is equivalent to the decomposition of signal variance in the frequency domain, whereas the bispectrum is the decomposition of signal skewness in the frequency domain. Therefore, the bispectrum can be used to describe the asymmetric and non-linear characteristics of signals.

Detection of Non-Linear Coupling Features in Signals
Because the bispectrum contains phase information, it can be used to determine whether non-linearity exists in the signal by detecting the quadratic phase coupling. Assume there is a sine wave x(n) = cos(k 0 n), with the period of N = 2π k 0 . The Fourier transform determined by the principal value sequence is: where k = ±1, X(k) = 0. The bispectrum is estimated as B x (k 1 , k 2 ) = X(k 1 )X(k 2 )X * (k 1 + k 2 ). Since k 1 = ±1, k 2 = ±1 and k 1 + k 2 = ±1 cannot hold at the same time, B x (k 1 , k 2 ) is equal to zero. Assume there is a sinusoidal signal with the harmonic y(n) = cos(k 0 n) + cos(2k 0 n), with a period of N = 2π k 0 . The Fourier transform determined by the principal value sequence is: where k = ±1 or k = ±2, Y(k) = 0. The bispectrum is estimated as B y (k 1 , . Obviously, B y (k 1 , k 2 ) has six non-zeros (k 1 = 1, k 2 = −2), Compared with the case without a harmonic signal, the bispectrum of y(n) is coupled at 6 points. It can be concluded that the frequency coupling features can be obtained by the bispectrum in the analysis of non-linear signals.

Immunity to Gaussian Noise
Assume there are two vibration signals, x(n) = e(n), y(n) = s(n) + e(n), where s(n) is the basic signal containing the second harmonic and e(n) is Gaussian noise. The traditional power-spectrum estimation method is used to analyze the two vibration signals, Obviously, the ability of power-spectrum estimation to detect the second harmonic is related to the signal-to-noise ratio (SNR). The lower the SNR, the lower the accuracy of the second harmonic will be. If bispectral analysis is used, i.e., In general, if the probability density function of noise is approximately symmetric about the longitudinal axis, B e (ω 1 , ω 2 ) can be ignored, so then the Equations (7) and (8) become: It can be seen that, even in the condition of low SNR, the SNR can be improved by bispectrum analysis, which is more conducive to the detection of non-linear signals. Zhou et al. [38] showed that the bispectrum analysis method could obtain more effective detection than the traditional Fourier transform method in the case of low SNR for signals mixed with random or Gaussian noise.

Calculation of the Bispectrum
The methods of high-order spectrum estimation can be divided into two categories: parametric and non-parametric. The advantages of non-parametric high-order spectrum estimation include ease of implementation and availability of the Fourier transformation (FFT). Moreover, non-parametric algorithms can be divided into direct and indirect methods [29]. For this paper, direct non-parametric bispectrum estimation is used.
Assuming that the random signal {x(n)} has N sample points, i.e., x(1), x(2), . . . , x(N), the implementation steps are as follows: (1) Sample segmentation: The N-point data is divided into K segments. Each segment contains M sample values, i.e., N = KM. M should meet the general length requirement of FFT. The data in the j th segment are recorded as x (j) (n), j = 1, 2, . . . , K. If needed, zeros can be added at the end of each segment to meet the length requirement. As with power spectrum estimation, the segments can overlap to meet the requirements of frequency resolution and estimation variance.
(2) Calculation of discrete Fourier coefficients: where, λ = 0, 1, . . . , M 2 − 1. (3) Calculation of the third order cumulant: where, ∆ 0 = f s N 0 . (4) Calculation of the bispectrum average: A contradiction between estimation variance and frequency resolution has frequently been encountered in traditional spectrum estimation methods. In practical applications, the final frequency resolution, f s N , is selected on the premise of balanced estimation variance and frequency resolution.

Numerical Simulation
In this section, the finite element (FE) method was used to study the non-linear dynamic characteristics of a beam with a breathing crack.

Finite Element (FE) Modeling of a Cantilever Beam with a Breathing Crack
As illustrated in Figure 1, a three-dimensional, rectangular section cantilever beam was built using FE software ANSYS with solid45 elements which are defined by eight nodes having three degrees of freedom at each node. The geometric parameters of the beam, including the length (L), width (B), height (H), and the material parameters, including elastic modulus (E), Poisson's ratio (ν), and density (ρ), are shown in Table 1. The cantilever beam contained a unilateral penetrating breathing crack, with the depth of a and the distance from the fixed end of x c . A harmonic excitation F sin 2πωt was applied on the upper surface near the free end of the beam, and a sensing point was arranged on the lower surface near the free end to trace and record the dynamic responses of the beam. In the model, there were a total of 5120 elements, and the elements near the crack are refined to precisely capture the dynamic behaviors of the crack. The depth and position of the crack were signified by the dimensionless parameters p = a/H and q = x c /L, respectively. where, (4) Calculation of the bispectrum average: A contradiction between estimation variance and frequency resolution has frequently been encountered in traditional spectrum estimation methods. In practical applications, the final frequency resolution, s f N , is selected on the premise of balanced estimation variance and frequency resolution.

Numerical Simulation
In this section, the finite element (FE) method was used to study the non-linear dynamic characteristics of a beam with a breathing crack.

Finite Element (FE) Modeling of a Cantilever Beam with a Breathing Crack
As illustrated in Figure 1, a three-dimensional, rectangular section cantilever beam was built using FE software ANSYS with solid45 elements which are defined by eight nodes having three degrees of freedom at each node. The geometric parameters of the beam, including the length (L), width (B), height (H), and the material parameters, including elastic modulus (E), Poisson's ratio (ν), and density (ρ), are shown in Table 1. The cantilever beam contained a unilateral penetrating breathing crack, with the depth of a and the distance from the fixed end of c x . A harmonic excitation sin 2 F t   was applied on the upper surface near the free end of the beam, and a sensing point was arranged on the lower surface near the free end to trace and record the dynamic responses of the beam. In the model, there were a total of 5120 elements, and the elements near the crack are refined to precisely capture the dynamic behaviors of the crack. The depth and position of the crack were signified by the dimensionless parameters / p a H  and / c q x L  , respectively.  For simplicity, an open crack model is often used to ignore the non-linear characteristics induced by cracks. However, the open crack model is not consistent with actual cases, and studies have shown that crack closure should be taken into account for more accurate representation of structural responses. Gudmundson et al. [7] pointed out that the change of structural natural frequency, due to open cracks, is much smaller than those due to  For simplicity, an open crack model is often used to ignore the non-linear characteristics induced by cracks. However, the open crack model is not consistent with actual cases, and studies have shown that crack closure should be taken into account for more accurate representation of structural responses. Gudmundson et al. [7] pointed out that the Sensors 2021, 21, 1177 6 of 20 change of structural natural frequency, due to open cracks, is much smaller than those due to breathing cracks. Moreover, the open-crack model may underestimate the severity of a fatigue crack, resulting in increased risks to safety. Qian et al. [8] found that the difference in the displacement responses between intact and cracked beams decreases when an open-close crack model is used.
In this work, determining the opening and closing behaviors of the breathing crack was considered a local contact problem. There are three common contact algorithms used by ANSYS, this time using the generating function algorithm. Taking one crack surface as the target surface and the other as the contact surface, the contact behaviors of cracks were studied using the surface-to-surface contact method, as illustrated in Figure 2. Specifically, the target surface was allowed to penetrate the contact surface, but the contact surface was not allowed to penetrate the target face. During the vibration of the cantilever beam, there were three states of the breathing crack: (a) the crack was completely open, and no nodes on the target and contact surfaces were in contact; (b) the crack was completely closed, and all nodes on the target and contact crack surfaces were in contact; (c) in the transition state featuring pressure release or pressure accumulation, the nodes on the target and contact crack surfaces were in partial contact.
Sensors 2020, 20, x FOR PEER REVIEW 6 of 21 breathing cracks. Moreover, the open-crack model may underestimate the severity of a fatigue crack, resulting in increased risks to safety. Qian et al. [8] found that the difference in the displacement responses between intact and cracked beams decreases when an openclose crack model is used. In this work, determining the opening and closing behaviors of the breathing crack was considered a local contact problem. There are three common contact algorithms used by ANSYS, this time using the generating function algorithm. Taking one crack surface as the target surface and the other as the contact surface, the contact behaviors of cracks were studied using the surface-to-surface contact method, as illustrated in Figure 2. Specifically, the target surface was allowed to penetrate the contact surface, but the contact surface was not allowed to penetrate the target face. During the vibration of the cantilever beam, there were three states of the breathing crack: (a) the crack was completely open, and no nodes on the target and contact surfaces were in contact; (b) the crack was completely closed, and all nodes on the target and contact crack surfaces were in contact; (c) in the transition state featuring pressure release or pressure accumulation, the nodes on the target and contact crack surfaces were in partial contact.

Non-Linear Dynamic Characteristics of the Cantilever Beam Containing a Breathing Crack
The dynamic characteristics of the cracked-beam structure, including modal frequencies and mode shapes, were obtained by modal analysis. Figure 3, below, shows the first six modal frequencies and mode shapes of the beam with p = 50% and q = 30%. It shows that the cracks in modes 1, 4, and 5 were in a closed state. In mode 4, in particular, the target surface penetrated the contact surface. The crack in modes 2, 3, and 6 was open, where modes 2 and 3 correspond to the transition state, and it is difficult to judge whether mode 6 was in transition or in a fully opened state.

Non-Linear Dynamic Characteristics of the Cantilever Beam Containing a Breathing Crack
The dynamic characteristics of the cracked-beam structure, including modal frequencies and mode shapes, were obtained by modal analysis. Figure 3, below, shows the first six modal frequencies and mode shapes of the beam with p = 50% and q = 30%. It shows that the cracks in modes 1, 4, and 5 were in a closed state. In mode 4, in particular, the target surface penetrated the contact surface. The crack in modes 2, 3, and 6 was open, where modes 2 and 3 correspond to the transition state, and it is difficult to judge whether mode 6 was in transition or in a fully opened state.
To diagnose damage in rotating machinery, the use of harmonics as typical signs of non-linearity of structural responses has been common, not only in research but also in engineering practice. However, structural non-linear characteristics such as harmonics caused by a breathing crack are difficult to identify using a linear damage-detection method.
Hence, bispectral analysis was employed to identify the relationships between the harmonic components in the signal, which could be used to identify the harmonic coupling in the response and to estimate the degree of non-linearity.
A harmonic excitation of ω = 1 2 f 1 , F = 500 N was applied at the position shown above in Figure 1. A dynamic time-history calculation of the structure was carried out, the step size was 5 × 10 −4 s, and the total analysis time was 2.548 s. The acceleration response at the sensing point was extracted and bispectral analysis performed.
The result of the damaged beam is shown below in Figure 4a; Figure 4b shows the bispectrum of the beam without damage. Figure 4b shows only the main peak A, caused by the excitation frequency coupling at (57 Hz, 57 Hz), and peaks D 1 and D 2 , excited by the load without coupling at (57 Hz, 0 Hz) and (0 Hz, 57 Hz).
Compared with Figure 4b, several harmonics can be seen in Figure 4a in terms of the coupling subpeaks of the excitation frequency and the second harmonic B 1 and B 2 at (57 Hz, 114 Hz) and (114 Hz, 57 Hz), and the self-coupling subpeak of the second harmonic C at (114 Hz, 114 Hz), along with other harmonic components.

Non-Linear Dynamic Characteristics of the Cantilever Beam Containing a Breathing Crack
The dynamic characteristics of the cracked-beam structure, including modal frequencies and mode shapes, were obtained by modal analysis. Figure 3, below, shows the first six modal frequencies and mode shapes of the beam with p = 50% and q = 30%. It shows that the cracks in modes 1, 4, and 5 were in a closed state. In mode 4, in particular, the target surface penetrated the contact surface. The crack in modes 2, 3, and 6 was open, where modes 2 and 3 correspond to the transition state, and it is difficult to judge whether mode 6 was in transition or in a fully opened state. To diagnose damage in rotating machinery, the use of harmonics as typical signs of non-linearity of structural responses has been common, not only in research but also in engineering practice. However, structural non-linear characteristics such as harmonics caused by a breathing crack are difficult to identify using a linear damage-detection method.
Hence, bispectral analysis was employed to identify the relationships between the harmonic components in the signal, which could be used to identify the harmonic coupling in the response and to estimate the degree of non-linearity.
A harmonic excitation of ω = ，F = 500 N was applied at the position shown above in Figure 1. A dynamic time-history calculation of the structure was carried out, the step size was 5 × 10 s, and the total analysis time was 2.548 s. The acceleration response at the sensing point was extracted and bispectral analysis performed.
The result of the damaged beam is shown below in Figure 4a; Figure 4b shows the bispectrum of the beam without damage. Figure 4b shows only the main peak A, caused by the excitation frequency coupling at (57Hz, 57 Hz), and peaks and , excited by the load without coupling at (57 Hz, 0 Hz) and (0 Hz, 57 Hz).  Compared with Figure 4b, several harmonics can be seen in Figure 4a in terms of the coupling subpeaks of the excitation frequency and the second harmonic and at (57 Hz, 114 Hz) and (114 Hz, 57 Hz), and the self-coupling subpeak of the second harmonic C at (114 Hz, 114 Hz), along with other harmonic components.
This shows that harmonic and frequency coupling can be used to identify non-linear damage and to provide damage-identification indicators for the development of a structural health diagnosis and crack identification methods. By using bispectral analysis, not only can the harmonics of each order, caused by non-linearity, be observed, but also the coupling relations among frequencies.
In addition to the acceleration response, the bispectrum can also analyze displacement, velocity, and stress responses. The corresponding results are shown below in Figure  5a-c, respectively. In the figures, the second-order harmonics can be observed, although their amplitudes exhibit clear deviations. This shows that harmonic and frequency coupling can be used to identify non-linear damage and to provide damage-identification indicators for the development of a structural health diagnosis and crack identification methods. By using bispectral analysis, not only can the harmonics of each order, caused by non-linearity, be observed, but also the coupling relations among frequencies.
In addition to the acceleration response, the bispectrum can also analyze displacement, velocity, and stress responses. The corresponding results are shown below in Figure 5a-c, respectively. In the figures, the second-order harmonics can be observed, although their amplitudes exhibit clear deviations.

The Influence of Crack Depth on Non-Linear Dynamic Features
Because crack propagation is unavoidable during structural service life, crack-depth identification is critical. For this reason, understanding the relationship between non-linear dynamic characteristics and crack-depth variation is crucial.
We assumed that the crack was located at q = 0.1, but the crack depth was set to p =

The Influence of Crack Depth on Non-Linear Dynamic Features
Because crack propagation is unavoidable during structural service life, crack-depth identification is critical. For this reason, understanding the relationship between non-linear dynamic characteristics and crack-depth variation is crucial.
We assumed that the crack was located at q = 0.1, but the crack depth was set to p = 0.1, 0.2, 0.3, 0.4, and 0.5. The excitation frequency was ω = 1 2 f 1 . As explained in [39], a 1/2 super harmonic resonance was expected to occur at this frequency, and here, frequency components were easily fully generated. The magnitude of the excitation was F = 500 N.
The acceleration responses were extracted for bispectral analysis, giving rise to threedimensional bispectral images, as shown in Figure 6. The amplitudes of the second harmonic at each crack depth were extracted, with the variation trends shown below in Figure 7. It was apparent that the amplitude of the second harmonic increased monotonically in proportion to the increase in crack depth. Therefore, the amplitude of the second harmonic in the bispectrum can be used as a damage-sensitive feature to indicate the degree of damage of the cracked beam.

The Influence of Crack Depth on Non-Linear Dynamic Features
Because crack propagation is unavoidable during structural service life, crack-depth identification is critical. For this reason, understanding the relationship between non-linear dynamic characteristics and crack-depth variation is crucial.
We assumed that the crack was located at q = 0.1, but the crack depth was set to p = 0.1, 0.2, 0.3, 0.4, and 0.5. The excitation frequency was . As explained in [39], a 1/2 super harmonic resonance was expected to occur at this frequency, and here, frequency components were easily fully generated. The magnitude of the excitation was F = 500 N. The acceleration responses were extracted for bispectral analysis, giving rise to threedimensional bispectral images, as shown in Figure 6. The amplitudes of the second harmonic at each crack depth were extracted, with the variation trends shown below in Figure 7. It was apparent that the amplitude of the second harmonic increased monotonically in proportion to the increase in crack depth. Therefore, the amplitude of the second harmonic in the bispectrum can be used as a damage-sensitive feature to indicate the degree of damage of the cracked beam.

The Influence of Crack Location on Non-Linear Dynamic Features
For the beam with a breathing crack, as studied for this paper, the local stress field in the beam changed with the position of the damage, which in turn changed the non-linear dynamic features of the structure. Assuming a fixed crack depth of p = 0.5, the crack location was set to q = 0.05, 0.1, 0.3, 0.5, and 0.7, respectively. The excitation frequency and magnitude were kept at

The Influence of Crack Location on Non-Linear Dynamic Features
For the beam with a breathing crack, as studied for this paper, the local stress field in the beam changed with the position of the damage, which in turn changed the non-linear dynamic features of the structure. Assuming a fixed crack depth of p = 0.5, the crack location was set to q = 0.05, 0.1, 0.3, 0.5, and 0.7, respectively. The excitation frequency and magnitude were kept at ω = 1 2 f 1 and F = 500 N, respectively. The three-dimensional bispectral images obtained are shown above in Figures 8a-e and 9, below, shows the variation trend of the second harmonic amplitude with the crack position.
This result indicates that the closer the crack location is to the fixed end, the more significant is the structural non-linearity.

The Influence of Harmonic Excitation Frequency on Non-Linear Dynamic Features
To better study the non-linearity caused by the periodic movement of the breathing crack, the excitation frequency was adjusted to be ω = 1 3 f 1 , 1 2 f 1 , f 1 , 2 f 1 , and 3 f 1 , to investigate the variation features of non-linear dynamic characteristics. The crack location and depth were kept at q = 0.3 and p = 0.5, respectively. Acceleration responses were extracted, and bispectral analysis was then performed to obtain the three-dimensional bispectrum images, shown below in Figure 10.
It can be seen that, when the structure was subjected to harmonic excitation, multiple harmonic components in the response could be analyzed by the bispectral method, including the coupling of excitation frequencies, n-order harmonic components, coupling of n-order harmonic components and m-order harmonic components.
The coupling of excitation frequencies and the second harmonic was obvious in each bispectral image. When the excitation frequency was an integer times (n = 1, 2, 3) or one over an integer times (n = 1/3, 1/2) the first-order frequency of the structure, the harmonic components of the cracked beam could be fully excited. Frequency coupling at ω = 1 3 f 1 , 1 2 f 1 , f 1 , and 2 f 1 , was significant. Specifically, the coupling effect at ω = 1 3 f 1 was more complex, whereas ω = f 1 may have induced too much uncertainty to the structural vibration. Therefore, ω = 1 2 f 1 and ω = 2 f 1 could be chosen as the optimal excitation frequencies.   To better study the non-linearity caused by the periodic movement of the breathing crack, the excitation frequency was adjusted to be f , 1 2 f , and 1 3 f , to investigate the variation features of non-linear dynamic characteristics. The crack location and depth were kept at q = 0.3 and p = 0.5, respectively. Acceleration responses were extracted, and bispectral analysis was then performed to obtain the three-dimensional bispectrum images, shown below in Figure 10.
It can be seen that, when the structure was subjected to harmonic excitation, multiple harmonic components in the response could be analyzed by the bispectral method, including the coupling of excitation frequencies, n-order harmonic components, coupling of n-order harmonic components and m-order harmonic components.
The coupling of excitation frequencies and the second harmonic was obvious in each bispectral image. When the excitation frequency was an integer times (n = 1, 2, 3) or one over an integer times (n = 1/3, 1/2) the first-order frequency of the structure, the harmonic components of the cracked beam could be fully excited. Frequency coupling at

The Influence of Harmonic Excitation Magnitude on Non-Linear Dynamic Features
Another important parameter likely to cause differences in bispectral analysis results was the magnitude of excitation. Assuming that the crack location and depth were q = 0.3 and p = 0.5, respectively, and the excitation frequency was ω = 1 2 f 1 , six excitation amplitudes, of F = 50 N, 100 N, 250 N, 400 N, 750 N and 1000 N, respectively, were studied.
The bispectrum analysis shows results similar to those above in Figure 10b, and the variation trend of the second harmonic amplitude with the harmonic excitation magnitude is shown below in Figure 11. In particular, in order to study the influence of harmonic exci-tation magnitude on non-linear dynamic features, the amplitudes of the second harmonic were not normalized.
Sensors 2020, 20, x FOR PEER REVIEW 13 of 21 Figure 11. Variation trend of the second harmonic amplitude with the harmonic excitation magnitude.
It was evident that when the excitation magnitude was magnified by 20 times, from 50 N to 1000 N, the non-linear characteristics of the bispectrum, the amplitude of the second harmonic also increased.

The Influence of Noise on Non-Linear Dynamic Characteristics
In actual measurement scenarios, noise from different environmental sources (such as thermal noise, mechanical relaxation, and electrical noise from the data-acquisition system [40]) inevitably affects the measured response data. Noise is likely to affect the representation of non-linear features, and can even drown out the non-linear features of the response. V.A. Krysko-Jr. concluded in [41] that with the increase in the external transverse load, even a relatively small intensity of white noise has an essential influence on the character of structural vibrations. Theoretically, the bispectrum analysis method used for this paper was immune to Gaussian noise. This section aims to verify the immunity of bispectrum analysis to the influence of noise.
If the crack location and depth are q = 0.3 and p = 0.5, respectively, the excitation frequency and magnitude would be and F = 500 N, respectively. Gaussian white noise was introduced into the acceleration response to simulate the effect of environmental noise. The noise intensity was characterized by the SNR, defined as 10lg s n P SNR P  (14) where s P represents the effective power of the signal and n P represents the effective power of the noise.
The acceleration responses were extracted and bispectral analysis was used to obtain three-dimensional bispectral images, as shown below in Figure 12. It can be clearly seen that when the SNR is greater than 30 dB, there is no significant change in the bispectrum. However, when the SNR reaches 10 dB, disturbances caused by Gaussian white noise begin to appear. The amplitudes of the second harmonic, under different noise levels, are shown below in Figure 13. It is clear that although the amplitude of the second harmonic of the bispectrum fluctuates around the noiseless signal due to the influence of noise, the characterization of non-linear dynamic features was not seriously impaired, even under the SNR of 10 dB. Such observation indicates that bispectral analysis has strong noise resistance in processing non-linear signals. It was evident that when the excitation magnitude was magnified by 20 times, from 50 N to 1000 N, the non-linear characteristics of the bispectrum, the amplitude of the second harmonic also increased.

The Influence of Noise on Non-Linear Dynamic Characteristics
In actual measurement scenarios, noise from different environmental sources (such as thermal noise, mechanical relaxation, and electrical noise from the data-acquisition system [40]) inevitably affects the measured response data. Noise is likely to affect the representation of non-linear features, and can even drown out the non-linear features of the response. V.A. Krysko-Jr. concluded in [41] that with the increase in the external transverse load, even a relatively small intensity of white noise has an essential influence on the character of structural vibrations. Theoretically, the bispectrum analysis method used for this paper was immune to Gaussian noise. This section aims to verify the immunity of bispectrum analysis to the influence of noise.
If the crack location and depth are q = 0.3 and p = 0.5, respectively, the excitation frequency and magnitude would be ω = 1 2 f 1 and F = 500 N, respectively. Gaussian white noise was introduced into the acceleration response to simulate the effect of environmental noise. The noise intensity was characterized by the SNR, defined as where P s represents the effective power of the signal and P n represents the effective power of the noise. The acceleration responses were extracted and bispectral analysis was used to obtain three-dimensional bispectral images, as shown below in Figure 12. It can be clearly seen that when the SNR is greater than 30 dB, there is no significant change in the bispectrum. However, when the SNR reaches 10 dB, disturbances caused by Gaussian white noise begin to appear. The amplitudes of the second harmonic, under different noise levels, are shown below in Figure 13. It is clear that although the amplitude of the second harmonic of the bispectrum fluctuates around the noiseless signal due to the influence of noise, the characterization of non-linear dynamic features was not seriously impaired, even under the SNR of 10 dB. Such observation indicates that bispectral analysis has strong noise resistance in processing non-linear signals.

Experimental Validation
In this section, two experiments were carried out to investigate the effects of crack and excitation parameters on the non-linear dynamic features of the cantilever beam with a breathing crack. Figure 13. Variation trend of the second harmonic amplitude with signal-to-noise ratio (SNR).

Experimental Validation
In this section, two experiments were carried out to investigate the effects of crack and excitation parameters on the non-linear dynamic features of the cantilever beam with a breathing crack.

Experiment 1
A cantilever beam with constant rectangular section, as shown below in Figure 14, was made with a length of L = 0.6 m, cross-section dimensions H × B = 0.019 m × 0.012 m, material density 7750 kg/m 3 , elastic modulus E = 1.84 × 10 11 Pa, and Poisson's ratio ν = 0.3. Figure 13. Variation trend of the second harmonic amplitude with signal-to-noise ratio (SNR).

Experimental Validation
In this section, two experiments were carried out to investigate the effects of crack and excitation parameters on the non-linear dynamic features of the cantilever beam with a breathing crack.

Experiment 1
A cantilever beam with constant rectangular section, as shown below in Figure 14  An accelerometer was mounted on the beam's surface 90 mm from the free end, and hammer excitation was applied 45 mm from the free end. The acceleration response was amplified by a charge amplifier and then collected by the UTel data-acquisition system at a sampling frequency of 12,800 Hz.
This experiment mainly investigated the influence of crack parameters on the nonlinear dynamic features of the cracked beam. The crack parameters are shown in Table 2. An accelerometer was mounted on the beam's surface 90 mm from the free end, and hammer excitation was applied 45 mm from the free end. The acceleration response was amplified by a charge amplifier and then collected by the UTel data-acquisition system at a sampling frequency of 12,800 Hz.
This experiment mainly investigated the influence of crack parameters on the non-linear dynamic features of the cracked beam. The crack parameters are shown in Table 2. The acceleration responses under different crack parameters were extracted, and the amplitude of the second harmonic components with different crack depths and positions were analyzed using the bispectrum. Figure 15a,b show the variation trend of the second harmonic amplitude with the crack depth and location, respectively. As with the numerical results, the amplitude of the second harmonic of the bispectrum increased along with the increase in crack depth. On the other hand, as the crack moved closer to the fixed end, the amplitude of the second harmonic of the bispectrum increased continuously.
amplitude of the second harmonic components with different crack depths and positions were analyzed using the bispectrum. Figure 15a,b show the variation trend of the second harmonic amplitude with the crack depth and location, respectively. As with the numerical results, the amplitude of the second harmonic of the bispectrum increased along with the increase in crack depth. On the other hand, as the crack moved closer to the fixed end, the amplitude of the second harmonic of the bispectrum increased continuously.

Experiment 2
A steel cantilever beam with the dimensions 380 mm × 20 mm × 10 mm was made, having the material elastic modulus of 193 GPa and Poisson's ratio of 0.29. The dynamic displacement response was measured in the experiment. To generate crack breathing in the beam, a v-shaped notch 1 mm in depth was first created on one side of the beam, 308 mm from the free end. Cyclic load was then applied to extend the crack depth to 5 mm.
Ranging from the fixed end to 4 mm from the free end, 11 measurement points were uniformly arranged along the centerline of the beam, as shown above in Figure 16. An excitation was applied 10 mm from the fixed end. Natural frequencies of different orders of the structure were captured within the ranges 0-2 kHz, 2-20 kHz, and 20-50 kHz, respectively. The first three orders of the natural frequencies measured 1 f = 70.63 Hz, 2 f = 688.70 Hz, and 3 f = 1787.19 Hz, respectively. .

Experiment 2
A steel cantilever beam with the dimensions 380 mm × 20 mm × 10 mm was made, having the material elastic modulus of 193 GPa and Poisson's ratio of 0.29. The dynamic displacement response was measured in the experiment. To generate crack breathing in the beam, a v-shaped notch 1 mm in depth was first created on one side of the beam, 308 mm from the free end. Cyclic load was then applied to extend the crack depth to 5 mm.
Ranging from the fixed end to 4 mm from the free end, 11 measurement points were uniformly arranged along the centerline of the beam, as shown above in Figure 16. An excitation was applied 10 mm from the fixed end. Natural frequencies of different orders of the structure were captured within the ranges 0-2 kHz, 2-20 kHz, and 20-50 kHz, respectively. The first three orders of the natural frequencies measured f 1 = 70.63 Hz, f 2 = 688.70 Hz, and f 3 = 1787.19 Hz, respectively.

Experiment 2
A steel cantilever beam with the dimensions 380 mm × 20 mm × 10 mm was made, having the material elastic modulus of 193 GPa and Poisson's ratio of 0.29. The dynamic displacement response was measured in the experiment. To generate crack breathing in the beam, a v-shaped notch 1 mm in depth was first created on one side of the beam, 308 mm from the free end. Cyclic load was then applied to extend the crack depth to 5 mm.
Ranging from the fixed end to 4 mm from the free end, 11 measurement points were uniformly arranged along the centerline of the beam, as shown above in Figure 16. An excitation was applied 10 mm from the fixed end. Natural frequencies of different orders of the structure were captured within the ranges 0-2 kHz, 2-20 kHz, and 20-50 kHz, respectively. The first three orders of the natural frequencies measured 1 f = 70.63 Hz, 2 f = 688.70 Hz, and 3 f = 1787.19 Hz, respectively. This experiment mainly investigated the influence of excitation parameters on the non-linear dynamic features of the cracked beam. The excitation parameters are shown above in Table 3. The acceleration responses under different excitation parameters were extracted, and the amplitude of the second harmonic varying with the excitation frequency and magnitude were subjected to bispectral analysis. A laser-scanning vibrometer, Polytec PSV-400, was used for non-contact measurement of the displacements. PZT piezoceramics were used to excite the structure, and the magnitude of the excitation was proportional to the output amplitude of the piezoceramics generator. Thus, changes in the magnitude of excitation were achieved by changing the voltage. Table 3. Harmonic excitation parameters.

Group
Excitation Frequency (Hz) Voltage (v) 1 1/2 0.125 This experiment mainly investigated the influence of excitation parameters on the non-linear dynamic features of the cracked beam. The excitation parameters are shown above in Table 3. The acceleration responses under different excitation parameters were extracted, and the amplitude of the second harmonic varying with the excitation frequency and magnitude were subjected to bispectral analysis. A laser-scanning vibrometer, Polytec PSV-400, was used for non-contact measurement of the displacements. PZT piezoceramics were used to excite the structure, and the magnitude of the excitation was proportional to the output amplitude of the piezoceramics generator. Thus, changes in the magnitude of excitation were achieved by changing the voltage.  Figure 17 below shows the amplitude of the second harmonic varying with voltage, which was directly correlated with the magnitude of excitation. The experimental results were similar to those of the numerical simulation. When the voltage was magnified by 48 times, from 0.125 V to 6 V, the amplitude of the second harmonic also increased.
quency and magnitude were subjected to bispectral analysis. A laser-scanning vibrometer, Polytec PSV-400, was used for non-contact measurement of the displacements. PZT piezoceramics were used to excite the structure, and the magnitude of the excitation was proportional to the output amplitude of the piezoceramics generator. Thus, changes in the magnitude of excitation were achieved by changing the voltage.  Figure 17 below shows the amplitude of the second harmonic varying with voltage, which was directly correlated with the magnitude of excitation. The experimental results were similar to those of the numerical simulation. When the voltage was magnified by 48 times, from 0.125 V to 6 V, the amplitude of the second harmonic also increased.  Figure 18 shows the bispectrum under different excitation frequencies. In the figure, the second harmonic of the bispectrum under different excitation frequencies can be observed. When the frequency is one over an integer times of the first-order frequency of the structure (e.g., n = 1/3, 1/2), the non-linearity can be characterized well and the coupling of other frequencies can be observed. However, when the harmonic excitation frequency is an integer times the first-order frequency (e.g., n = 1, 2, 3), the second harmonic is not obvious, especially under the frequency of 3 , where the frequency coupling is more  Figure 18 shows the bispectrum under different excitation frequencies. In the figure, the second harmonic of the bispectrum under different excitation frequencies can be observed. When the frequency is one over an integer times of the first-order frequency of the structure (e.g., n = 1/3, 1/2), the non-linearity can be characterized well and the coupling of other frequencies can be observed. However, when the harmonic excitation frequency is an integer times the first-order frequency (e.g., n = 1, 2, 3), the second harmonic is not obvious, especially under the frequency of 3 f 1 , where the frequency coupling is more complex. Therefore, ω = 1/2 f 1 is expected to generate both obvious and explicit frequency coupling phenomena. 0.637%, respectively. When the excitation frequency was 1 1 2 f , the error was smaller, and when the excitation frequency was 1 1 3 f , the error was slightly larger.
However, the magnitude of error did not affect the non-linear characterization of the bispectrum. The immunity of bispectrum analysis to noise influence is demonstrated well. Figure 19. Variation trend of the second harmonic amplitude with noise.

Conclusions
For this paper, the non-linear dynamic characteristics of a cantilever beam containing a breathing crack were studied using bispectral analysis. Numerical simulation and experimental verification show that, compared with traditional damage identification techniques, the bispectral analysis method demonstrates sensitivity to the non-linear features of the cracked beam and can effectively overcome the limitations of existing linear damage detection methods in non-linear damage detection. Several conclusions can be drawn: (1) The bispectral analysis method can not only characterize the non-linear features of crack structures, but can also reveal the coupling relationships among frequencies.
(2) The non-linear feature of the bispectrum, in terms of the amplitude of the second harmonic, is proportional to crack depth, so the amplitude of the second harmonic of the bispectrum can be used as a damage-sensitive feature to quantify the severity of cracks. With the crack position effect, the closer the crack is to the fixed end, the larger is the amplitude of the bispectral second harmonic, and the more obvious are the non-linear features.
(3) As demonstrated in the numerical simulation and experimental results, the harmonic components in the dynamic response signal are the most abundant, and the non-linear characteristic effect is best when the excitation frequency  is  However, the magnitude of error did not affect the non-linear characterization of the bispectrum. The immunity of bispectrum analysis to noise influence is demonstrated well.

Conclusions
For this paper, the non-linear dynamic characteristics of a cantilever beam containing a breathing crack were studied using bispectral analysis. Numerical simulation and experimental verification show that, compared with traditional damage identification techniques, the bispectral analysis method demonstrates sensitivity to the non-linear features of the cracked beam and can effectively overcome the limitations of existing linear damage detection methods in non-linear damage detection. Several conclusions can be drawn: (1) The bispectral analysis method can not only characterize the non-linear features of crack structures, but can also reveal the coupling relationships among frequencies.
(2) The non-linear feature of the bispectrum, in terms of the amplitude of the second harmonic, is proportional to crack depth, so the amplitude of the second harmonic of the bispectrum can be used as a damage-sensitive feature to quantify the severity of cracks. With the crack position effect, the closer the crack is to the fixed end, the larger is the amplitude of the bispectral second harmonic, and the more obvious are the non-linear features. (3) As demonstrated in the numerical simulation and experimental results, the harmonic components in the dynamic response signal are the most abundant, and the non-linear characteristic effect is best when the excitation frequency ω is 1 2 f 1 . The amplitude of the second harmonic increases with the increase of the excitation. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data sharing not applicable No new data were created or analyzed in this study. Data sharing is not applicable to this article.