Evaluation of fatigue cracks using nonlinearities of acousto-ultrasonic waves acquired by an active sensor network

Without loss of generality, consider a one-dimensional rod guiding a longitudinal wave mode. By ignoring the wave attenuation and dispersion, the wave can be defined by

$$\begin{eqnarray} &&\rho \frac{{\partial }^{2}u(x,t)}{\partial {t}^{2}}=\frac{\partial \sigma }{\partial x}, \end{eqnarray} \tag{ 1 }$$

where ρ is the medium density, u(x,t) the particle displacement at location x and moment t (abbreviated as u in what follows; u is namely a captured AU wave signal), and σ the normal stress. Taking into account the medium nonlinearity, the nonlinear Hooke's law reads [29]

$$\begin{eqnarray} &&\sigma =E \varepsilon (1+\beta \varepsilon +\cdots ), \end{eqnarray} \tag{ 2 }$$

where E is the modulus of elasticity of the medium, ε the normal strain and β the second-order nonlinear elastic coefficient (β is called acoustic nonlinearity factor hereinafter). Substituting equations (2) into (1) and neglecting the terms higher than the second order (due to their minute contribution to σ), it gives rise to

$$\begin{eqnarray} &&\rho \frac{{\partial }^{2}u}{\partial {t}^{2}}=E\frac{{\partial }^{2}u}{\partial {x}^{2}}+E \beta \frac{\partial u}{\partial x}\frac{{\partial }^{2}u}{\partial {x}^{2}}. \end{eqnarray} \tag{ 3 }$$

To solve equation (3), a perturbation theory [30] is applied, in which u is re-defined as

$$\begin{eqnarray} &&u={u}_{0}+{u}_{1}, \end{eqnarray} \tag{ 4 }$$

where u₀ represents the incident AU wave, a zero-order term governing the linear components of u, and u₁ is the first-order perturbation accounting for the nonlinear components of u. Substituting equations (4) into (3) leads to

$$\begin{eqnarray} &&\fl \rho \frac{{\partial }^{2}{u}_{0}}{\partial {t}^{2}}-E\frac{{\partial }^{2}{u}_{0}}{\partial {x}^{2}}+\rho \frac{{\partial }^{2}{u}_{1}}{\partial {t}^{2}}-E\frac{{\partial }^{2}{u}_{1}}{\partial {x}^{2}}\nonumber\\ &&-~E \beta \frac{\partial ({u}_{0}+{u}_{1})}{\partial x}\frac{{\partial }^{2}({u}_{0}+{u}_{1})}{\partial {x}^{2}}=0. \end{eqnarray} \tag{ 5 }$$

Because the nonlinear term in equation (5) (i.e., $\rho \frac{{\partial }^{2}{u}_{1}}{\partial {t}^{2}}-E\frac{{\partial }^{2}{u}_{1}}{\partial {x}^{2}}-E \beta \frac{\partial ({u}_{0}+{u}_{1})}{\partial x}\frac{{\partial }^{2}({u}_{0}+{u}_{1})}{\partial {x}^{2}}$) is sufficiently small when compared to the linear term (i.e., $\rho \frac{{\partial }^{2}{u}_{0}}{\partial {t}^{2}}-E\frac{{\partial }^{2}{u}_{0}}{\partial {x}^{2}}$), one has,

$$\begin{eqnarray} &&\rho \frac{{\partial }^{2}{u}_{0}}{\partial {t}^{2}}=E\frac{{\partial }^{2}{u}_{0}}{\partial {x}^{2}}. \end{eqnarray} \tag{ 6 }$$

By substituting equations (6) into (5), one obtains

$$\begin{eqnarray} &&\rho \frac{{\partial }^{2}{u}_{1}}{\partial {t}^{2}}=E\frac{{\partial }^{2}{u}_{1}}{\partial {x}^{2}}+E \beta \frac{\partial ({u}_{0}+{u}_{1})}{\partial x}\frac{{\partial }^{2}({u}_{0}+{u}_{1})}{\partial {x}^{2}}. \end{eqnarray} \tag{ 7 }$$

Provided u₁ << u₀ (nonlinear attributes in the medium are considerably weaker than linear attributes), equation (7) reduces to

$$\begin{eqnarray} &&\rho \frac{{\partial }^{2}{u}_{1}}{\partial {t}^{2}}=E\frac{{\partial }^{2}{u}_{1}}{\partial {x}^{2}}+E \beta \frac{\partial {u}_{0}}{\partial x}\frac{{\partial }^{2}{u}_{0}}{\partial {x}^{2}}. \end{eqnarray} \tag{ 8 }$$

Assume u₀ to be a harmonic wave with a mono-frequency of ω, as

$$\begin{eqnarray} &&{u}_{0}={\Gamma }_{1}\cos (k x-\omega t), \end{eqnarray} \tag{ 9 }$$

where k is the wavenumber, and Γ₁ represents the magnitude of the first-order harmonic wave, which is called the primary wave mode in what follows. With equation (9), the corresponding perturbation solution to equation (8) can be obtained as

$$\begin{eqnarray} &&{u}_{1}=\frac{\beta }{8}{\Gamma }_{1}^{2}{k}^{2}x\cos (2 k x-2 \omega t), \end{eqnarray} \tag{ 10a }$$

or, if denoting $\frac{\beta }{8}{\Gamma }_{1}^{2}{k}^{2}x$ as Γ₂,

$$\begin{eqnarray} &&{u}_{1}={\Gamma }_{2}\cos (2 k x-2 \omega t). \end{eqnarray} \tag{ 10b }$$

Equation (10) reveals that the nonlinear component u₁ with a magnitude of Γ₂, characterizing the nonlinear attributes of the captured AU wave signal, has a nominal frequency of 2ω—double that of the excitation frequency. This implies that the medium nonlinearity distorts an incident mono-frequency AU wave, resulting in an energy shift to higher frequencies (integer multiples of the excitation frequency), a phenomenon defined as acoustic nonlinearity [18]. The wave components at higher frequencies are collectively referred to as higher-order harmonics of the primary wave mode, and in particular the one with the double driving frequency (2ω) is the second harmonic wave mode. From equation (10), one has

$$\begin{eqnarray} &&\beta =\frac{8}{{k}^{2}x}\frac{{\Gamma }_{2}}{{\Gamma }_{1}^{2}}. \end{eqnarray} \tag{ 11 }$$

Therefore, β, the acoustic nonlinearity factor hereinafter, links the magnitudes of the primary (Γ₁) and second harmonic wave modes (Γ₂) to the propagation distance (x).

The acoustic nonlinearity factor defined by equation (11) was derived for longitudinal waves propagating in a one-dimensional rod. Regulated by a scaling parameter (γ), such a factor can be extended to Lamb waves in plates as detailed elsewhere [20, 29, 30], as

$$\begin{eqnarray} &&\beta =\frac{8}{{k}^{2}x}\frac{{\Gamma }_{2}}{{\Gamma }_{1}^{2}}\gamma . \end{eqnarray} \tag{ 12 }$$

Because γ depends in a complex manner on the material and geometric parameters of the medium, a relative acoustic nonlinearity factor, β', is introduced in this study, in recognition of the fact that a medium has an unchanged γ at a given propagation distance of AU waves before and after the occurrence of fatigue damage. This factor reads

$$\begin{eqnarray} &&{\beta }^{\prime}=\frac{{\Gamma }_{2}}{{\Gamma }_{1}^{2}}. \end{eqnarray} \tag{ 13 }$$

This relative acoustic nonlinearity factor introduces additional benefits to fatigue damage characterization: (i) in principle, the characterization is independent of a benchmark structure and a baseline signal, because the process is based on the detection of any abnormal increase in β' due to fatigue damage, which is considerably minute in a pristine benchmark (without fatigue damage). In other words, the nonlinearity manifested in a captured AU wave signal contributed by the material itself is insignificant compared with that arising from the fatigue damage (to be demonstrated in section 3.3), and therefore once any singular increase in β' detected, fatigue damage is deemed existent; (ii) it minimizes the influence of the material and geometric diversity of the medium on AU nonlinearity; and (iii) it avoids the involvement of wave propagation distance as seen in equation (12).

To assure considerable generation and accumulation of the desired second harmonic of AU waves, the primary and second harmonic wave modes should conform to the following cardinal prerequisites: (i) synchronism: both the phase and the group velocities of the primary wave mode match those of the second harmonic wave mode, respectively, thus guaranteeing the shift of AU energy from the primary to the second harmonic modes over propagation distance (thus the second harmonic mode becomes accumulative) [19, 31]; (ii) non-zero power flux: the primary mode is of the same type as the second harmonic mode (e.g., both are symmetric or both are anti-symmetric), thus ensuring a non-zero power transfer from the primary to the second harmonic modes. Any two wave modes meeting the above two conditions are referred to as a synchronous pair in this study; in contrast, power transfer from the primary to the second harmonic modes could not be efficiently achieved if the above two conditions were not satisfied.

Figure 1 shows the dispersion curves of AU Lamb waves in an aluminum plate, calculated using DISPERSE^®. The mode pair (S₁,S₂) at (f = 3.57 MHz⋅mm, 2f = 7.14 MHz⋅mm)—the first-order (denoted by S₁) and second-order (denoted by S₂) symmetric Lamb modes, as highlighted in the figure—meets the prerequisites and forms a synchronous pair. In addition, at 3.57 MHz⋅mm and 7.14 MHz⋅mm, the S₁ and S₂ modes propagate at the highest speeds among all available modes, respectively, thus standing out from the other modes and avoiding spurious contributions from other wave modes to the second harmonic accumulation. Such a merit excludes other candidate synchronous pairs such as (S₂,S₄). For the selected synchronous pair, the magnitude of S₁ at 3.57 MHz⋅mm corresponds to Γ₁ in equation (13), while the magnitude of S₂ at 7.14 MHz⋅mm corresponds to Γ₂.

Zoom In Zoom Out Reset image size

The aluminum plate underwent a fatigue processing on a digitally controlled MTS testing platform (858 Mini Bionix), subject to a sinusoidal tensile load of 5 Hz. It took about 200 000 cycles to produce a barely visible fatigue crack (BVFC), circa 4 mm in length from the tip of the initiator, as seen in figure 3. Upon generation of the BVFC, ten circular PZT wafers (nominal diameter: 5 mm each) were surface-mounted on the fatigued plate, to form five actuator–sensor paths A_i − S_i (i = I,II,III,IV,V), each having different od (0, 10, 20, 30, 40 mm, respectively). For each path, the distance between the actuator and the sensor was kept at 200 mm. Each path was instrumented with a signal generation and acquisition system [32]. 16-cycle Hanning-window-modulated sinusoidal tone bursts at a central frequency of 800 kHz were generated in MATLAB and downloaded to an arbitrary waveform generation unit (Agilent E1441). The signal was amplified with a high-frequency power amplifier (CIPRIAN US-TXP-3) to 80 V_p−p and applied on each actuator of the five sensing paths to generate AU waves. With the current plate thickness (4.5 mm), the excitation frequency of 800 kHz ensured generation of S₁ as the primary wave mode (0.8 MHz × 4.5 mm = 3.57 MHz⋅mm), which constituted a synchronous pair with S₂. The response signals (including S₁ and S₂) were acquired by the PZT sensor of each sensing path through a signal digitizer (Agilent E1438) at a sampling rate of 40 MHz.

Zoom In Zoom Out Reset image size

As an example, the time domain signal acquired by path A_I − S_I is exhibited in figure 4(a). Often the captured signals present highly multimodal and dispersive natures, from which no straightforward signal features associated with fatigue damage can possibly be extracted directly. It is noteworthy that a tiny BVFC does not incur noticeable changes in linear macroscopic wave scattering, and thus no additional wave packets scattered from the BVFC or mode conversion can be detected in the signals. This once again accentuates that identification relying on linear macroscopic wave scattering may lose its effectiveness when used to detect fatigue damage. The signal in figure 4(a) was subsequently processed with the short-time Fourier transform (STFT) to obtain its time–frequency spectrogram, shown in figure 4(b), highlighting two major energy concentrations at the excitation frequency (circa 800 kHz) and twice the frequency (circa 1.6 MHz), respectively. The primary and second harmonic wave modes were extracted at their respective frequencies in figure 4(b) and re-constructed to the time domain, combined in figure 5 for comparison. From figure 5, Γ₁ of S₁ and Γ₂ of S₂ were determined to calculate β' using equation (13). Note that Γ₂ is only 5% of Γ₁, indicating the weakness of the second harmonic mode compared with the primary mode from which it is generated and accumulated. It is noteworthy that a certain discrepancy in the arrival time of the maximum energy between S₁ and S₂ can be noticed in figure 5, which can be attributed to possible errors while experimentally capturing two wave modes (as explained S₂ is insignificant compared with S₁) and extracting signal features using the designated signal processing tool.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Repeating the above process, values of β' for all five sensing paths were calculated, contributing to a correlation curve between the relative acoustic nonlinearity factor of a particular sensing path and its od from the BVFC, after normalization with the one calculated via A_I − S_I, as shown in figure 6 (the solid line). In the figure, different values of od were normalized with regard to the AU wavelength (od/λ), making it possible to extend the discussion to general circumstances with other excitation frequencies.

Zoom In Zoom Out Reset image size

Dynamic FE analysis was conducted in parallel using the commercial FE package ABAQUS^®/EXPLICIT. The same aluminum plate was meshed at full scale using three-dimensional eight-node brick elements. To be consistent with the experiment configuration, encastre boundary conditions were applied on the left and right plate edges. A surface crack 4 mm long from the tip of the initiator and 2.25 mm deep (half of the intact plate thickness), parallel to the left/right plate edge, was modelled. In the vicinity of the crack, the mesh was particularly refined with an element size less than 0.25 mm, thus guaranteeing at least 15 elements were allocated per wavelength of the primary wave mode. Eight elements were allocated through the plate thickness.

A contacting-pair interaction was imposed on the two surfaces of the crack, which had a gap of 0.2 mm in between, connecting the two surfaces in a contact manner and enforcing the equality of their degrees of freedom. Such a restriction allowed the breathing behaviour of a fatigue crack when the primary wave mode traversed it, properly simulating CAN. A three-dimensional nonlinear constitutive relationship involving the third-order elastic (TOE) tensor was introduced to the material property definition, so as to model the intrinsic material nonlinearity and the geometric nonlinearity [33] as well, through a user subroutine VUMAT in ABAQUS^®, as fully reported in the authors' other work [29]. A pre-developed piezoelectric actuator/sensor model [32] was recalled to model the PZT wafers, in which uniform in-plane radial displacement constraints were applied on the FE nodes of the actuator model to generate the primary wave mode, which is consistent with the motion pattern of the S₁ mode, which predominantly has in-plane displacements. The same excitation condition as that used in the experiment was applied. AU wave signals were acquired by calculating the strains at the places where sensors were positioned.

To facilitate comparison with experimental results, the combined primary and second harmonic wave modes extracted from the signal captured by the same sensing path A_I − S_I are shown in figure 7, from which β' was ascertained using equation (13). Compared with the corresponding experimental results shown in figure 5, excellent coincidence between two wave modes in terms of the arrival time of the maximum energy is observed, as a result of the absence of any measurement errors in simulation.

Zoom In Zoom Out Reset image size

With all five sensing paths, the dotted line in figure 6 presents the relationship between β' of a particular sensing path and its od/λ from the crack obtained from simulation. Quantitative agreement is noticed between the results from experiment and from the enhanced FE analysis using TOE for simulating nonlinear AU waves.

The correlation β' − od/λ in figure 6 reveals that the relative acoustic nonlinearity factor defined by equation (13) captured via A_I − S_I (i.e., the path having the smallest offset distance from the fatigue crack) presents the highest value; with the increase in od, β' decreases monotonously and reaches its minimum when od/λ is approximately 6 (corresponding to the case that od is circa 30 mm); in particular, a remarkable decrease in β' is observed within od/λ < 4 (od: 20 mm), whereas no further significant decrease can be noted when od/λ > 4.

The above correlation was established when the length of the sensing path (between A_i and S_i) was 200 mm. To further scrutinize this correlation at different distances of AU wave propagation, modelling and simulations were repeated for different sensing path lengths, ranging from 100 to 400 mm with an increment of 100 mm. The accordingly obtained β' − od/λ curves are compared in figure 8, where no substantial variations among the results obtained by sensing paths of different lengths can be observed, provided their od from the crack remained the same. The largest β' can always be captured when the sensing path has the shortest od, and it decreases gradually with increasing od and subdues once od/λ reaches approximately 4. This affirms the hypothesis made in section 2.2 that compared with the accumulative material nonlinearity as wave propagation, the one incurred due to fatigue damage dominates the overall nonlinearity of an AU wave signal (i.e., the nonlinearity of an AU wave signal contributed by the material itself is insignificant compared with that arising from the fatigue damage). The correlation in figure 6 rationalizes the algorithm developed in this study for quantitative characterization of fatigue damage.

Zoom In Zoom Out Reset image size

In DDI, the inspection area of the structure is two-dimensionally meshed first using a dense set of spatial points, and projected to a two-dimensional pixelated image (or map) with each pixel corresponding exclusively to a spatial point. The value of each image pixel (called field value hereinafter) is in some way linked to the presence probability of fatigue damage at the spatial point correlated with this pixel. A fatigue damage zone can thus be graphically defined by highlighting those pixels with remarkably high field values, providing a quantitative and detailed depiction about the damage (e.g., size and orientation).

Considering an active PZT sensor network with N_Total PZT elements to cover the inspection area, each PZT in the network can serve as either the actuator or the senor in terms of the dual piezoelectric effects. Given a sensing path in the network PZT_i − PZT_j(i,j = 1,2,⋯N_Total,i ≠ j), with the coordinates of actuator PZT_i and sensor PZT_j being (x_i,y_i) and (x_j,y_j), respectively, the field value at pixel (x,y) perceived by PZT_i − PZT_j (denoted by I(x,y)_i−j) is defined, based on the conclusion drawn from figure 6, as

$$\begin{eqnarray} &&I(x,y)_{i-j}={\beta }_{{}_{i-j}}^{\prime}\left [\frac{\varsigma -{R}_{i j}(x,y)}{\varsigma -1}\right ], \end{eqnarray} \tag{ 14 }$$

where ${\beta }_{{}_{i-j}}^{\prime}$ is the relative acoustic nonlinearity factor acquired via sensing path PZT_i − PZT_j, ς is a scaling parameter controlling the size of the effective distribution area, and R_ij(x,y) is a weight to regulate the area of influence from the fatigue damage on a particular sensing path in the sensor network as detailed elsewhere [34, 38], which reads

$$\begin{eqnarray} &&\fl {R}_{i j}(x,y)\nonumber\\ =~\left \{ \begin{array}{@{}ll@{}} \displaystyle \frac{\sqrt{(x-{x}_{i})^{2}+(y-{y}_{i})^{2}}+\sqrt{(x-{x}_{j})^{2}+(y-{y}_{j})^{2}}}{\sqrt{({x}_{i}-{x}_{j})^{2}+({y}_{i}-{y}_{j})^{2}}}&\displaystyle \mathrm{when}\hspace{0.167em} \lt \varsigma \\ \displaystyle \varsigma &\displaystyle \mathrm{when}\hspace{0.167em} \geq \varsigma . \end{array}\right. \normalsize \end{eqnarray} \tag{ 15 }$$

With equation (14), the field value at each pixel can be constructed by PZT_i − PZT_j, leading to a probability image. In the case that the inspection area is free of damage, all the field values are low (practically it is not zero due to noise interference). In contrast, with fatigue damage, the sensing paths can be affected by the nonlinearity due to fatigue damage, more or less subject to the od between the fatigue damage and this sensing path, and the field values at those pixels contained in the fatigue damage zone are elevated pronouncedly.

It is important to note that a sensing path perceives the damage near it only, an attribute due to the quick descent of β' as with an increase in od, as described in figures 6 and 8. This high sensitivity creates the feasibility of identifying multiple fatigue damage zones, because such a trait of β' makes the field value defined by equation (14) highly inert to distant damage (to be demonstrated in section 5).

In accordance with equation (14), each sensing path in the sensor network creates a probability image (called source image) [39], each reflecting the prior perceptions as to the fatigue damage in the structure from the perspective of the sensing path that creates such an image. Each source image contains, however, not only features related to damage but unwanted components such as measurement noise and uncertainties. Aimed at screening unwanted noise while strengthening a fatigue damage-incurred singularity, an image fusion scheme based on the arithmetic mean was introduced, amalgamating available source images at the pixel level and constructing a superimposed image (ultimate resulting image). It is a collective consensus as to fatigue damage from the entire sensor network. The scheme is defined as

$$\begin{eqnarray} &&I(x,y)\biggr \vert _{\mathrm{ultimate}}=\frac{1}{N}\sum _{{i=1,j=1\atop (i\not = j)}}^{{N}_{\mathrm{total}}}I(x,y)_{i-j}, \end{eqnarray} \tag{ 16 }$$

where I(x,y)|_ultimate is the field value at pixel (x,y) in the ultimate resulting image. Upon image fusion, the damage-incurred singularity (commonality in individual source images) stands out and noise (random information in source images) is filtered. Identified fatigue damage zone can thus be illuminated at image pixels with higher field values.

An aluminum plate (measuring 380 × 400 × 4.5 mm³) with four through-thickness rivet holes for bolt connection (Φ10 mm each; see figure 9 and table 1 for geometric parameters) was fatigued using a digitally controlled MTS testing platform (858 Mini Bionix) subject to a sinusoidal tensile load of 5 Hz. To facilitate initiation of micro-fatigue cracks and expedite their deterioration, a small notch with a sharp tip was inscribed at the edge of two rivet holes randomly chosen (#1 and #2), respectively. It took about 500 000 cycles to produce two barely visible fatigue cracks, as photographed in figure 10: about 5 mm long initiated from Rivet Hole #1, and circa 3 mm from #2. Eight PZT wafers (Φ: 5 mm, thickness: 0.5 mm each), denoted by T_i(i = 1,2,...,8) with respective coordinates tabulated in table 1, were surface-mounted on the fatigued plate to form an active PZT sensor network, as shown schematically in figure 9. This sensor network offered 7 × 8 = 56 sensing paths, which were instrumented with the signal generation and acquisition system introduced in section 3.1. The diagnostic signal (Hanning-window-modulated 16-cycle sinusoid tone bursts at 800 kHz, conforming to the mode selection criterion in section 2.3) was generated and the response AU wave signals were captured with the system.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With captured signals, β' was calculated using equation (13) for each sensing path, and a source image was subsequently constructed using equation (14). All source images were fused in terms of equation (16), leading to the ultimate resulting image, as displayed in figure 11. In the image, the presence probability of fatigue damage is quantitatively calibrated in terms of the greyscale, where the darker the greyscale at a pixel the greater the presence possibility of fatigue damage at the spatial node represented by that pixel. For easy comparison, the two real fatigue cracks are integrated in the image. In the figure, two regions with higher greyscale (i.e., higher field values) can be observed, covering Rivet Hole #1 and #2, as well as their vicinities, implying a higher probability for the existence of fatigue cracks therein, although artefacts due to measurement noise can also be seen in the image which perturb the evaluation results from the reality.

Zoom In Zoom Out Reset image size

It is noteworthy that the identified results in the image correspond to the two fatigue cracks initiated from rivet holes, rather than the rivet holes themselves or the fatigue crack initiators at the hole edges, although the highlighted areas with higher greyscale are considerably greater than the actual sizes of the two fatigue cracks. In fact, the two highlighted zones describe the vicinity of the two fatigue cracks. That is because, upon fatigue testing, the vicinities of two fatigue cracks became plastically deformed, which elevates β' measured by the sensing paths nearby. As explained in section 2.2, the proposed approach explores the abnormal increase in β' due to a fatigue crack only, and it is thus independent of a benchmark structure, in which the intrinsic material nonlinearity, geometric nonlinearity and gross damage (e.g., a rivet hole in this study) would not incur abnormal changes in β'. Thus, once the increase in β' is detected, it is attributable to the existence of fatigue damage.

Furthermore, a threshold, κ, was applied to strengthen the identified fatigue damage by screening measurement noise. κ is a pre-set percentage of the maximum field value of the ultimate resulting image. Any field value which is less than such a threshold is forced to be zero, and otherwise retained (a smaller threshold means more information is included in the ultimate resulting image). Shown in figure 12 are the resulting images when different thresholds (κ = 0.5,0.6,0.7,0.8) are applied. With the increase in κ, artefacts from noise are considerably suppressed and the fatigue damage zones become obvious. In particular, with κ = 0.8 and ς = 1.01, two fatigue damage zones are manifested clearly in figure 12(d).

Zoom In Zoom Out Reset image size