# Nonlinear acoustic modeling and measurements during the fatigue process in metals

Wenhan Lyu, Xianmei Wu, Wei-Jiang Xu

<figure><img src="image_1.png" /></figure>

## To cite this version:

Wenhan Lyu, Xianmei Wu, Wei-Jiang Xu. Nonlinear acoustic modeling and measurements during the fatigue process in metals. Materials, 2019, 12 (4), pp.607. 10.3390/ma12040607. hal-03152515

HAL Id: hal-03152515

https://hal.science/hal-03152515

Submitted on 25 Feb 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Article

# Nonlinear Acoustic Modeling and Measurements
during the Fatigue Process in Metals

Wenhan Lyu<sup>1,2,3</sup>, Xianmei Wu<sup>1,2,*</sup> and Weijiang Xu<sup>3</sup>

Received: 11 January 2019; Accepted: 15 February 2019; Published: 18 February 2019

<figure><img src="image_3.png" /></figure>

check for
updates

**Abstract:** The nonlinear spring model combined with dislocation dipole theory was applied to describe the acoustic nonlinearity during the fatigue process in metals. The spring stiffness changes with fatigue degree. For the early stage, spring stiffness approaches infinity, and the heavier nonlinearity mainly results from the increase of dislocation density. Further fatigue leads to the occurrence of micro-cracks, during which spring stiffness begins to decrease. Abundant micro-crack sprouting accelerates the crack's expansion, and spring stiffness drops quickly, which causes the obvious decline in the transmitted harmonic amplitudes. Solutions obtained from the nonlinear wave equation with dislocation terms were added into the spring model. Varying spring stiffness was chosen for simulating the fatigue process. Then, nonlinear harmonic variation during this process was observed, which was classified into three stages: (I) the early dislocation fatigue stage; (II) the micro-crack sprouting stage; (III) the crack expansion stage. Nonlinear acoustic measurements were carried out on an aluminum alloy specimen during its fatigue process until cracks could be seen clearly. Harmonic variations in experiments can also be classified into the same three stages as the numerical results, which provides a theoretical and experimental reference for fatigue evaluation in metals using the nonlinear acoustic method.

**Keywords:** nonlinear spring model; dislocation dipoles; fatigue evaluation

## 1. Introduction

Metal materials are widely used in practical engineering applications, but long-term in-service work makes them prone to failure due to fatigue damages. The linear acoustic method is invalid for detecting fatigue damages of materials, while the finite amplitude sound wave can produce obvious waveform distortion when it propagates through a medium with fatigue damages, which has been proven to be an effective method for evaluating fatigue damages in metal materials [1–5].

Nonlinear acoustic theory in solids established by researchers such as Landau, Murnaghan, and Goldberg provides the basis for nonlinear ultrasonic testing [6,7]. Based on this, theoretical models were proposed to explain the mechanism of the nonlinear acoustic effect caused by micro-damages. For early stage fatigue damages are usually described by the dislocation string [8] or dislocation dipole models [9,10], while these are not suitable for heavier damages owing to the sprouting of micro-cracks. The PM (proposed by Preisach and Mayergoyz) model [11] and contact interface model [12–15] were proposed to explain the nonlinear acoustic mechanism caused by micro-cracks or imperfect interfaces, through which nonlinear effects such as harmonics, difference frequency waves, resonance frequency drift, and slow dynamics can be simulated. Most of the research works have been based on the

<footnote n="1">¹ State Key Laboratory of Acoustics, Institute of Acoustics, Chinese Academy of Sciences, Beijing 100190, China; lvwenhan@mail.ioa.ac.cn</footnote>

<footnote n="2">² University of Chinese Academy of Sciences, Beijing 100049, China</footnote>

<footnote n="3">³ Université Polytechnique Hauts-de-France, CNRS, Univ. Lille, YNCREA, Centrale Lille, UMR 8520IEMN-DOAE, F-59313 Valenciennes CEDEX 9, France; wei-jiang.xu@uphf.fr</footnote>

<footnote>* Correspondence: wuxm@mail.ioa.ac.cn; Tel.: +86-010-8254-7811</footnote>
above theoretical models. For example, reflection and transmission theories at an open or closed interface were established to study the contact nonlinearity of acoustics based on the contact interface model [16,17]. Some experiments [9,10,18,19] were carried out to observe the nonlinear properties of material delamination and microstructure changes caused by fatigue damages, which established the connection between micro-damages and acoustical nonlinearity. Nonlinear Lamb waves have been commonly applied to evaluate the fatigue damages or detect the bonding quality of solid boards [20,21]. Furthermore, some researchers [22] proved the effectiveness of the nonlinear acoustic method in estimating the strength degradation caused by initial damage of concrete materials. High temperature also leads to micro-damages, which can be measured through the ultrasonic nonlinear parameters to evaluate the microstructure changes during the high-temperature degradation process [23]. The spring model is usually used to study the characteristics of nonlinear acoustic waves at the solid rough contact interface to evaluate the interface boundary conditions [14].

Nonlinear acoustics are used in many research works for the non-destructive evaluation of metal materials, but the above research works are mainly focused on one of the stages during the fatigue process, such as early stage fatigue evaluation or crack detection. Actually, fatigue is a continuously changing process, so modeling the whole process is essential for its evaluation. In our work, the nonlinear spring model combined with dislocation dipoles theory was established for modeling the continuous fatigue process from the early dislocation stage to the crack expansion stage. Solutions obtained from the nonlinear wave equation with dislocations [24–27] were added into the nonlinear spring model. For the early stage, spring stiffness approaches infinity. Fatigue damage is mainly caused by dislocations in this stage. For the late stage damages, varying stiffness was chosen for modeling the different fatigue degree with the sprouting of micro-cracks. Once the crack begins to expand, spring stiffness declines quickly, and the amplitude of each order harmonic drops obviously. The nonlinear acoustic experimental system was built, and measurements were carried out during the whole fatigue process of the sample, which certified the numerical models.

## 2. Theoretical Models

### 2.1. Nonlinearity Due to Lattice Dislocations at Early Stage Fatigue

Early stage fatigue damage of metals mainly originates from dislocations. As the degree of fatigue further deepens, the dislocation density increases progressively, and the nonlinear acoustic effects becomes more obvious. The dislocation can be described by the dipole model [10], as shown in Figure 1.

<figure><img src="image_4.png" /><figcaption>Figure 1. The dislocation dipole model: (a) CT specimen (b) dislocation dipole.</figcaption></figure>

A compact tension (CT) specimen (Figure 1a) was chosen for modeling. Stress during the fatigue loading mainly concentrates on the initial grooving. At the early stage of the loading process, the appearance of dislocations is considered as pairs of dipoles, as sketched in Figure 1b: $(x, h)$ is the
position of one dislocation P in the pair relative to another one O at (0,0). The two dislocation lines are parallel and lie perpendicular to the "xy" plane. The unit length glide force [9,10] that acts between the two dislocations in the x direction is given as follows:

$$F_x = -\frac{Gb^2}{2\pi(1-v)} \frac{x(x^2-h^2)}{(x^2+h^2)^2} \quad (1)$$

where G is the shear modulus, v is Poisson's ratio, b is the magnitude of Burger's vector, and h represents the spacing between the two glide planes. Propagation of an ultrasonic wave in the material produces the strains of two distinguished origins: the elastic strain εl due to the lattice structure of the material and the strain εd associated with the motion of dislocation dipoles.

$$\varepsilon = \varepsilon_l + \varepsilon_d \quad (2)$$

The relationship between the elastic strain εl and stress σ is given by [9,28]:

$$\varepsilon_l = \frac{1}{A_2^l} \sigma - \frac{1}{2} \frac{A_3^l}{(A_2^l)^3} \sigma^2 + \dots \quad (3)$$

where A_2^l, A_3^l are the second- and third-order elastic constants, respectively, while the contribution of strain εd to the stress σ can be obtained from the dislocation motion equation [10]. Only the relative motion of the dipoles along the direction parallel to the glide planes is considered:

$$m_{eff} \frac{\partial^2 x}{\partial t^2} = F_x + b\sigma R \quad (4)$$

where m_{eff} is the effective mass of the dislocation dipoles, and R is the Schmid factor. Let ζ = x - h, which is the relative displacement of dislocation dipoles to the position h. Combining Equations (1) and (4) and expanding Equation (1) in the series method with ζ, the following is obtained:

$$m_{eff} \frac{\partial^2 \zeta}{\partial t^2} + \frac{Gb^2}{4\pi h^2(1-v)} \zeta - \frac{Gb^2}{8\pi h^3(1-v)} \zeta^2 + \dots = b\sigma R \quad (5)$$

Solving this equation and using the relationship εd = ΩΛdbζ, where, Ω is the conversion factor and Λd is the dislocation density, the following relationship is obtained:

$$\varepsilon_d = \frac{1}{A_2^d} \sigma - \frac{1}{2} \frac{A_3^d}{(A_2^d)^3} \sigma^2 + \dots \quad (6a)$$

with:

$$A_2^d = -\frac{G}{4\pi\Omega R \Lambda_d h^2 (1-v)}, \quad A_3^d = \frac{G}{4\pi\Omega^2 R \Lambda_d^2 h^3 (1-v)b} \quad (6b)$$

Substituting Equations (3) and (6a,b) into Equation (2) yields:

$$\varepsilon = \left( \frac{1}{A_2^l} + \frac{1}{A_2^d} \right) \sigma - \frac{1}{2} \left[ \frac{A_3^d}{(A_2^d)^3} + \frac{A_3^l}{(A_2^l)^3} \right] \sigma^2 + \dots \quad (7)$$

The inverse expression can be written in the following form approximately:

$$\sigma = A_2^l \varepsilon - \frac{1}{2} \left[ A_3^l + \frac{A_3^d (A_2^l)^3}{(A_2^d)^3} \right] \varepsilon^2 + \dots \quad (8)$$
Substituting Equation (8) and Equation (6b) into the following wave motion equation:

$$ \rho \frac{\partial^2 u}{\partial t^2} = \frac{\partial \sigma}{\partial x}, \quad \varepsilon = \frac{\partial u}{\partial x} \qquad (9) $$

we have:

$$ \frac{\partial^2 u}{\partial t^2} - c^2 \frac{\partial^2 u}{\partial x^2} = c^2 (\beta_l + \beta_d) \frac{\partial^2 u}{\partial x^2} \frac{\partial u}{\partial x} + \dots \qquad (10a) $$

with:

$$ \beta_l = - \frac{A_3^l}{A_2^l}, \quad \beta_d = - \frac{A_3^d (A_2^l)^2}{(A_2^d)^3} = \frac{16\pi^2 \Omega R^2 \Lambda_d h^3 (1 - v)^2 (A_2^l)^2}{G^2 b} \qquad (10b) $$

where $u$ is the total displacement, $\beta_l$, $\beta_d$ are the nonlinear coefficients due to the material lattices and dislocations, respectively, $c = \sqrt{A_2^l/\rho}$ is the wave velocity, and $\rho$ is the mass density of the material. Therefore, Equation (10a,b) is the nonlinear acoustic wave equation in materials with dislocations, from which the relationship between dislocation density and nonlinear harmonics can be established. The perturbation method is usually employed to solve this equation.

Assuming the following solution form that is expanded with the Mach number:

$$ u = \sum_{i=1}^{n} M^i u_i. \qquad (11) $$

where the Mach number $M = V_0/c$ and $V_0$ is the velocity amplitude of the vibration. Substituting Equation (11) into Equation (10a):

$$ \frac{1}{c^2} \sum_{i=1}^{n} \frac{\partial^2}{\partial t^2} (M^i u_i) = \sum_{i=1}^{n} \frac{\partial^2}{\partial x^2} (M^i u_i) + \beta \sum_{i=1}^{n} \frac{\partial}{\partial x} (M^i u_i) \frac{\partial^2}{\partial x^2} (M^i u_i) + \dots \qquad (12a) $$

with:

$$ \beta = \beta_l + \beta_d \qquad (12b) $$

Factoring the above equation in terms of $M^i$,

$$ M^1 : \frac{1}{c^2} \frac{\partial^2 u_1}{\partial t^2} = \frac{\partial^2 u_1}{\partial x^2}, \qquad (13a) $$

$$ M^2 : \frac{1}{c^2} \frac{\partial^2 u_2}{\partial t^2} = \frac{\partial^2 u_2}{\partial x^2} + \beta \frac{\partial u_1}{\partial x} \frac{\partial^2 u_1}{\partial x^2}, \qquad (13b) $$

$$ M^{\dots} : \dots \qquad (13c) $$

from which the fundamental and higher harmonics of the propagating wave can be resolved by the perturbation method.

The fundamental solution of Equation (13a) in the positive $x$ direction has the form:

$$ u_1 = A_0 \cos(kx - \omega t), \qquad (14) $$

where $A_0$ is the wave amplitude at the initial position (usually a driving source), $k$ is the wave number, and $\omega$ is the angular frequency.

For the second harmonic, substituting Equation (14) into Equation (13b),

$$ \frac{1}{c^2} \frac{\partial^2 u_2}{\partial t^2} - \frac{\partial^2 u_2}{\partial x^2} = \frac{1}{2} \beta A_0^2 k^3 \sin(2kx - 2\omega t) \qquad (15) $$
whose solution can be found to be,

$$u_2 = \frac{1}{8} \beta A_0^2 k^2 x \cos(2kx - 2\omega t) \quad (16)$$

In the same way, the former order solutions are used to determine the next order one:

$$u_3 = \beta^2 \left[ -\frac{1}{32} A_0^3 k^4 x^2 \cos(kx - \omega t) - \frac{3}{32} A_0^3 k^3 x \sin(kx - \omega t) + \frac{1}{32} A_0^3 k^4 x^2 \cos(3kx - 3\omega t) + \frac{5}{96} A_0^3 k^3 x \sin(3kx - 3\omega t) \right] \quad (17)$$

$$u_4 = \beta^3 \left[ \left( -\frac{1}{96} A_0^4 k^6 x^3 + \frac{1}{12} A_0^4 k^4 x \right) \cos(2kx - 2\omega t) - \frac{33}{512} A_0^4 k^5 x^2 \sin(2kx - 2\omega t) + \left( \frac{1}{96} A_0^4 k^6 x^3 - \frac{11}{512} A_0^4 k^4 x \right) \cos(4kx - 4\omega t) + \frac{5}{128} A_0^4 k^5 x^2 \sin(4kx - 4\omega t) \right] \quad (18)$$

and so on, to obtain $u_5, u_6, \dots, u_n$.

Substituting all the solutions into Equation (11) and combining the terms having the same harmonic frequency, we have,

$$u = \sum_n A_n \sin \left( nkx - n\omega t + \arctan \frac{c_n}{s_n} \right) \quad (19)$$

where $A_n = \sqrt{s_n^2 + c_n^2}$ and $\arctan \frac{c_n}{s_n}$ (for $n = 1, 2, 3, 4, \dots$) are the amplitude and phase of the nth harmonic wave component, respectively. $A_n, s_n, c_n$ are functions of propagation distance $x$ and wave amplitude $A_0$ at the initial position. The first four orders are discussed in our work and are given as follows:

$$s_1 = -\frac{3}{32} M^3 \beta^2 A_0^3 k^3 x, \quad c_1 = A_0 M - \frac{1}{32} M^3 \beta^2 A_0^3 k^4 x^2 \quad (20a)$$

$$s_2 = -M^4 \beta^3 \frac{33}{512} A_0^4 k^5 x^2, \quad c_2 = \frac{1}{8} M^2 \beta A_0^2 k^2 x + M^4 \beta^3 \left( -\frac{1}{96} A_0^4 k^6 x^3 + \frac{1}{12} A_0^4 k^4 x \right) \quad (20b)$$

$$s_3 = \frac{5}{96} M^3 \beta^2 A_0^3 k^3 x, \quad c_3 = \frac{1}{32} M^3 \beta^2 A_0^3 k^4 x^2 \quad (20c)$$

$$s_4 = \frac{5}{128} M^4 \beta^3 A_0^4 k^5 x^2, \quad c_4 = M^4 \beta^3 \left( \frac{1}{96} A_0^4 k^6 x^3 - \frac{11}{512} A_0^4 k^4 x \right) \quad (20d)$$

$$M = \frac{V_0}{c} = \frac{\omega A_0}{c} \quad (20e)$$

Through the above formulations, the relationship between the dislocation density and high order harmonics can be deduced.

## 2.2. Nonlinearity Due to Micro-Cracks at Late Stage Fatigue

While the fatigue of the metals is gradually strengthened, a further increase of the dislocation density leads to the occurrence of micro-cracks, which means a weak interface is produced in the material. The expansion of the micro-cracks will change the nonlinear behavior of the transmission wave, in both its transmitted energy and harmonics distribution. The acoustic nonlinearity mechanism in this case differs from that of the dislocation model. It can be established with the following nonlinear spring model, as shown in Figure 2.
<figure><img src="image_5.png" /><figcaption>Figure 2. Nonlinear spring model.</figcaption></figure>

A and B are two planes from either side of a crack interface. They are connected by the interfacial force, behaving like springs. Supposing that the displacement of an incident acoustic wave at plane A is $f(x - ct)$, and $g(x + ct)$, $h(x - ct)$ represents that of the reflected and transmitted wave, respectively. P is used to simulate the static force of the closed crack, and G is the joint dynamic force caused by the motion of the interface. The midpoint of the interface is taken as the original point "x₀", where $x_{0-}$, $x_{0+}$ represent the upper side A and the lower side B of the interface. $u(x_{0-}, t)$, $u(x_{0+}, t)$ are the displacement of sides A and B, respectively. The relationship between the interface spacing $y(t) = u(x_{0+}, t) - u(x_{0-}, t)$ and $G(y)$, $F(t)$ is given as the following [16,17],

$$ \rho c \dot{y}(t) = 2G(y) - 2F(t) \qquad (21) $$

where $F(t)$ is caused by the incident ultrasonic wave,

$$ F(t) = \left( -\rho c^2 \frac{\partial f(x - ct)}{\partial x} \right)_{x=x_0} \qquad (22) $$

From Equation (21), the variation of the crack pitch $y(t)$ is related to both the joint dynamic force and ultrasonic force, as expressed in Equation (22). $G(y)$ is usually nonlinear, whose relationship with interface spacing is given by the nonlinear spring model.

$$ G(y) = -K_1 y + K_2 y^2 + \dots \qquad (23) $$

where $K_1$, $K_2$ are the linear and nonlinear stiffness of the spring, respectively, which are related to the interface conditions. Considering the former dislocation model, the incident wave was already nonlinear because of the dislocation damage during the early stage before it propagated through the micro-cracks' area, so it can be obtained from the former results at position $x_0$ in the following form according to Equation (19):

$$ f(x - ct) = \sum_n A_n \sin(nkx_0 - n\omega t + \varphi_n) \qquad (24a) $$

with:

$$ \varphi_n = \arctan \frac{c_n}{s_n}, n = 1, 2, 3, \dots \qquad (24b) $$

Substitute Equations (23) and (24a,b) into Equation (21), then,

$$ \dot{y} + \frac{2K_1}{\rho c} y - \frac{2K_2}{\rho c} y^2 - \dots = 2 \sum_n n\omega A_n \cos(nkx_0 - n\omega t + \varphi_n) \qquad (25) $$

The solution of the above equation is composed of two parts: one linear part $y_1$ and a series of nonlinear parts $y_2, \dots$

$$ y = y_1 + y_2 + \dots, \quad y_1 \gg y_2 \gg \dots \qquad (26) $$
Then,

$$ \dot{y}_1 + \dot{y}_2 + \frac{2K_1}{\rho c}(y_1 + y_2 + \dots) - \frac{2K_2}{\rho c}(y_1 + y_2 + \dots)^2 - \dots = 2 \sum_n n\omega A_n \cos(nkx_0 - n\omega t + \varphi_n) \quad (27) $$

The linear part $y_1$ is easy to obtain by solving Equation (28):

$$ \dot{y}_1 + \frac{2K_1}{\rho c} y_1 = 2 \sum_n n \omega A_n \cos(nkx_0 - n\omega t + \varphi_n) \quad (28) $$

Then,

$$ y_1 = 2 \sum_n \frac{n A_n}{\sqrt{4K_1^2/n^2 \omega^2 \rho^2 c^2 + n^2}} \cos(nkx_0 - n\omega t + \varphi_n + \varphi_n) \quad (29) $$

where $\varphi_n = \arctan \frac{n\omega\rho c}{2K_1}$, and $y_2$ can be resolved similarly from the following equation, while the small terms $y_1y_2, y_2^2, \dots$ are omitted.

$$ \dot{y}_2 + \frac{2K_1}{\rho c} y_2 = \frac{2K_2}{\rho c} y_1^2 \quad (30) $$

The harmonics of four orders are considered in the result,

$$ y = \frac{2A_1 \cos(kx_0 - \omega t + \varphi_1 + \varphi_1)}{\sqrt{4K_1^2/\omega^2 \rho^2 c^2 + 1}} + \frac{4A_2 \cos(2kx_0 - 2\omega t + \varphi_2 + \varphi_2)}{\sqrt{4K_1^2/\omega^2 \rho^2 c^2 + 4}} + \frac{6A_3 \cos(3kx_0 - 3\omega t + \varphi_3 + \varphi_3)}{\sqrt{4K_1^2/9\omega^2 \rho^2 c^2 + 9}} + \frac{8A_4 \cos(4kx_0 - 4\omega t + \varphi_4 + \varphi_4)}{\sqrt{4K_1^2/16\omega^2 \rho^2 c^2 + 16}} - \frac{2A_1^2 K_2 \sin(2kx_0 - 2\omega t + 2\varphi_1 + 2\varphi_2 + \varphi_2)}{(4K_1^2/\omega^2 \rho^2 c^2 + 1)\sqrt{4K_1^2 + 4\omega^2 \rho^2 c^2}} - \frac{8K_2 A_2^2 \sin(4kx_0 - 4\omega t + \varphi_4 + 2\varphi_2 + \varphi_4)}{(4K_1^2/\omega^2 \rho^2 c^2 + 4)\sqrt{4K_1^2 + 16\omega^2 \rho^2 c^2}} + \frac{4K_2 A_2 A_1}{\sqrt{(4K_1^2/\omega^2 \rho^2 c^2 + 1)(4K_1^2/\omega^2 \rho^2 c^2 + 4)}} \left( \frac{\cos(kx_0 - \omega t + \varphi_3 - \varphi_1 + \varphi_2)}{\sqrt{4K_1^2 + \omega^2 \rho^2 c^2}} + \frac{\cos(3kx_0 - 3\omega t + \varphi_3 + \varphi_1 + \varphi_1 + \varphi_2 + \varphi_3)}{\sqrt{4K_1^2 + 9\omega^2 \rho^2 c^2}} \right) + \dots \quad (31) $$

The transmitted wave through the spring interface can be obtained with the following equations [16,17],

$$ h(x - ct) = f(x - ct) + \frac{1}{2} y \left( t - \frac{x}{c} \right) \quad (32) $$

Then,

$$ h(x - ct) = A_1 \sin(kx_0 - \omega t + kx + \varphi_1) + A_2 \sin(2kx_0 - 2\omega t + 2kx + \varphi_2) + A_3 \sin(3kx_0 - 3\omega t + 3kx + \varphi_3) + A_4 \sin(4kx_0 - 4\omega t + 4kx + \varphi_4) + \frac{A_1 \cos(kx_0 - \omega t + kx + \varphi_1 + \varphi_1)}{\sqrt{4K_1^2/\omega^2 \rho^2 c^2 + 1}} + \frac{2A_2 \cos(2kx_0 - 2\omega t + 2kx + \varphi_2 + \varphi_2)}{\sqrt{4K_1^2/9\omega^2 \rho^2 c^2 + 9}} + \frac{3A_3 \cos(3kx_0 - 3\omega t + 3kx + \varphi_3 + \varphi_3)}{\sqrt{4K_1^2/\omega^2 \rho^2 c^2 + 4}} + \frac{4A_4 \cos(4kx_0 - 4\omega t + 4kx + \varphi_4 + \varphi_4)}{\sqrt{4K_1^2/16\omega^2 \rho^2 c^2 + 16}} - \frac{2A_1^2 K_2 \sin(2kx_0 - 2\omega t + 2kx + 2\varphi_1 + 2\varphi_2 + \varphi_2)}{(4K_1^2/\omega^2 \rho^2 c^2 + 1)\sqrt{4K_1^2 + 4\omega^2 \rho^2 c^2}} - \frac{4K_2 A_2^2 \sin(4kx_0 - 4\omega t + 4kx + 2\varphi_2 + 2\varphi_2 + \varphi_2)}{(4K_1^2/\omega^2 \rho^2 c^2 + 4)\sqrt{4K_1^2 + 16\omega^2 \rho^2 c^2}} + \frac{2K_2 A_2 A_1}{\sqrt{(4K_1^2/\omega^2 \rho^2 c^2 + 1)(4K_1^2/\omega^2 \rho^2 c^2 + 4)}} \left( \frac{\cos(kx_0 - \omega t + kx + \varphi_3 - \varphi_1 + \varphi_2)}{\sqrt{4K_1^2 + \omega^2 \rho^2 c^2}} + \frac{\cos(3kx_0 - 3\omega t + 3kx + \varphi_3 + \varphi_1 + \varphi_1 + \varphi_2 + \varphi_3)}{\sqrt{4K_1^2 + 9\omega^2 \rho^2 c^2}} \right) + \dots \quad (33) $$

The transmitted coefficient of each order harmonic can be obtained with the above solution,

$$ T_n = \sqrt{\frac{4K_1^2/n^2 \omega^2 \rho^2 c^2 + 2n^2 + 2n \sin \phi_n \sqrt{4K_1^2/n^2 \omega^2 \rho^2 c^2 + n^2}}{4K_1^2/n^2 \omega^2 \rho^2 c^2 + n^2}}, \quad n = 1, 2, 3, 4 \dots \quad (34) $$

## 2.3. Computation Results

The transmission coefficients are related to the stiffness $K_1$, as shown in Figure 3, the second-order terms are small enough to be ignored compared with the linear parts. The computation parameters of the aluminum alloy are chosen as an example, and the values of relevant parameters are referenced from the work of John H. Cantrell [9], which are given in Table 1.
<table><caption>**Table 1.** Computational parameters of the aluminum alloy.</caption><thead><tr><td>Parameters</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>A<sub>2</sub><sup>l</sup></td><td>109</td><td>GPa</td></tr><tr><td>A<sub>3</sub><sup>l</sup></td><td>-510</td><td>GPa</td></tr><tr><td>G</td><td>28.6</td><td>GPa</td></tr><tr><td>b</td><td>0.4</td><td>mm</td></tr><tr><td>Ω</td><td>0.33</td><td>/</td></tr><tr><td>R</td><td>0.33</td><td>/</td></tr><tr><td>h</td><td>5</td><td>nm</td></tr><tr><td>v</td><td>0.33</td><td>/</td></tr><tr><td>c</td><td>6250</td><td>m/s</td></tr><tr><td>ρ</td><td>2700</td><td>kg/m<sup>3</sup></td></tr><tr><td>f</td><td>1</td><td>MHz</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Relationship between the transmitted coefficient of each order harmonic (the fundamental wave T₁, the second-order harmonics T₂, the third-order harmonics T₃, the fourth-order harmonics T₄) and spring stiffness.</figcaption></figure>

It is obvious that the transmitted coefficient of each order harmonic approaches one, while $K_1 \to \infty$, usually $K_1 \ge 1 \times 10^{15} N/m^3$, which can be considered as full transmission, as shown in Region I: the early dislocation fatigue stage in Figure 3. Therefore, during this stage, acoustic nonlinearity mainly results from dislocations and wave energy transmitted completely without the effects of the weak interface, as described in the dislocation dipole model. Once the micro-cracks sprout, from $K_1 = 1 \times 10^{15} N/m^3$, spring stiffness begins to decline, and the energy of the transmitted wave drops gradually, as shown in Region II: micro-crack sprouting stage. A large amount of micro-crack sprouting leads to the rapid expansion of the crack length, which reduces the transmitted wave energy clearly; it is the last stage of the fatigue process in Region III: crack expansion stage. During the micro-crack sprouting stage, nonlinearity increased, owing to the increasing dislocation density being very weak compared to the effects of spring the stiffness decline with continuous fatigue loading. Relative results are given in Figure 4. Computation parameters are also given in Table 1. The initial amplitude of each order harmonic is chosen from the above dislocation model results at $x_0 = 50\text{mm}$, driving source amplitude $A_0 = 1\mu\text{m}$. Assuming that micro-cracks sprout while the dislocation density is about $3 \times 10^{15} m^{-3}$, which is the early fatigue stage, the harmonics variations are given in Figure 4a,c. After this stage, the spring stiffness in the computation model was taken from $\infty \to 0$, which means $1/K_1$ from $0 \to \infty$, as shown in Figure 4b,d; the dislocation density was assumed as constant during this process.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Amplitude variation of (a) the fundamental wave $A_1$, (b) the second-order harmonics $A_2$, (c) the third-order harmonics $A_3$, and (d) the fourth-order harmonics $A_4$ during the fatigue process, which is expressed by dislocation density for early stage and spring stiffness for the later stage.</figcaption></figure>

The above results show that acoustic nonlinearity is very weak for materials without damages, as given in Figure 4c, while the dislocation density is zero. However, with the increase of dislocation density, which means the deepening of fatigue damage, acoustic nonlinearity becomes much heavier. The fundamental wave amplitude is going to decrease, while harmonic amplitudes begin to rise due to the energy distribution law that energy transfers from the lower order signals to the higher orders, as shown in Region I in Figure 4a,c. Then, acoustic nonlinearity indicates the fatigue degree in materials during the early fatigue stage. Micro-cracks sprout with further fatigue loading, and the value of spring stiffness begins to drop, whose effects on each order harmonic are much heavier than that of the increasing dislocation density, so the harmonic amplitudes are going to decrease, as shown in Region II in Figure 4b,d. The energy of each order transmitted wave declines very quickly with the decreasing stiffness, owing to the rapid expansion of cracks in Region III.

### 3. Experiment Results and Discussion

Based on the above analysis for the theoretical models during the fatigue process, the following experiments were designed and carried out. The experimental setup employed is represented in Figure 5a. It was composed of: a signal generator (Tektronix AFG1022, Tektronix, Shanghai, China), a power amplifier (Model 2100L RF Power Amplifier, ENI, Rochester, NY, USA), two transducers (central frequency of 1.1 MHz for the emitting one and 2.3 MHz for the receiving one, respectively), and a digital oscilloscope (Agilent Technologies MSO7032A, Agilent Technologies Inc., Santa Clara, CA, USA).
The receiving transducer sensitivity at 1.1 MHz was about −11.4 dB, for a signal of 2.2 MHz (the second-order signal), 3.3 MHz (the third-order one), and 4.4 MHz (the fourth-order one), and the sensitivities were −2 dB, −6 dB, and −12 dB, respectively. The gains of the receiving transducers were compensated in the following high-order harmonic amplitude measurements.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) The experiment system. (b) The stress distribution during the fatigue test.</figcaption></figure>

Before the preparation of the fatigue specimen, the fatigue model of the aluminum alloy specimen was established using finite element software (COMSOL Multiphysics 5.3, COMSOL Inc., Stockholm, Sweden). As in the fatigue test, the same size model was established with the same material. The lower surface boundary was fixed, and the upper surface boundary was loaded with the maximum stress during the fatigue test. The total stress distribution of this specimen in the fatigue process was analyzed, and the results are given in Figure 5b. It was shown that the stress in the test piece was concentrated along the initial grooving. If the fatigue degree increased constantly, tiny cracks would start to sprout and gradually extend forward. According to the calculation results of the model, the fatigue test was carried out on the sample by selecting appropriate parameters, the process of which is shown in Figure 6a.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** (a) The fatigue test. (b) The fatigue sample with a crack.</figcaption></figure>

The fatigue test was carried out at room temperature using the stress-loaded control mode, during which the average fatigue loading force was 5 kN and the loading frequency was 5 Hz. After a certain period of loading, nonlinear acoustic measurements were carried out and then redetected after another period of loading, until micro-cracks sprouted in the sample. The loading cycle was 500, 1500, 3000, 5000, 10,000, 15,000, 20,000, 25,000, and 29,800, respectively. Before the fatigue loading, acoustic measurement was performed to get the result of the non-damaged sample as a comparison, which is given in Figure 7. While the loading cycle was about 5000 during the fatigue test, acoustic nonlinearity increased very clearly in the received signals, as shown in Figure 8. Micro-cracks were observed from
the end of the preset grooving after 20,000 cycles of loading, the results of which are given in Figure 9. Once micro-cracks sprouted, continuous loading accelerated their expansion. Fatigue loading was performed until the crack could be seen clearly in the sample, as shown in Figure 6b, whose loading cycle was about 29,800, and the measurement results are given in Figure 10.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a) Transmitted signals received from the fatigue sample (Loading Cycle 0) in the time-domain and (b) in the frequency-domain (excitation voltage 200 V).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** (a) Transmitted signals received from the fatigue sample (Loading Cycle 5000) in the time-domain and (b) in the frequency-domain (excitation voltage 200 V).</figcaption></figure>

Comparing the results in Figures 7a and 8a, the signal received from the sample after the fatigue test (Loading Cycle 5000) deformed much more heavily than that of the non-fatigue situation (Loading Cycle 0), which is much more obvious in the spectrum results in Figures 7b and 8b. After fatigue loading, which means the increase of dislocation density, as explained in the dislocation dipole model, harmonics caused by the damage increased also; therefore, the spectrum was much richer in harmonic components, and each-order harmonic could be seen more clearly than the result of the non-damage situation.

The difference is not obvious in the time-domain result in Figure 9a, except that the transmitted energy is a little lower compared with the result in Figure 8a. However, in the spectrum results in Figure 9b, not only the harmonic amplitudes are a little lower, but another frequency component whose frequency is between the fundamental wave and second-order harmonics appears, compared with the results in Figure 8b, which may be caused by the wave scattering of the crack tip or the spatial variation in the stiffness of the springs. These results are much more apparent in Figure 10. As micro-cracks expanded very quickly with continuous loading, the energy of the transmitted wave in this stage was much lower than that in the former situations, as given in the time-domain result, which
was consistent with the spring model result. In addition, more spectrum components appeared in the spectrum result, as marked with rectangles in Figure 10b. Because the crack surface was rough and curved, this may have led to more scattering waves and generated more other intermediate frequency components in the result. It is very complicated to describe in theoretical models and will be studied in our next work.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** (a) Transmitted signals received from the fatigue sample (Loading Cycle 20,000) in the time-domain and (b) in the frequency-domain (excitation voltage 200 V).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** (a) Transmitted signals received from the fatigue sample (Loading Cycle 29800) in the time-domain and (b) in the frequency-domain (excitation voltage 200 V).</figcaption></figure>

The above results are some typical measurements for explaining the acoustic nonlinearity variation during different fatigue stages. Detailed amplitude variation of each order harmonic with fatigue degree, which is described by the fatigue loading cycle, is shown in Figure 11. During each measurement, excitation voltages were taken as 150 V (Figure 11a,b), 175 V (Figure 11c,d), and 200 V (Figure 11e,f), respectively.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Amplitude variation of (a) the fundamental wave A₁, (b) the high order harmonics Aₙ (excitation voltage 150 V), (c) the fundamental wave A₁, (d) the high order harmonics Aₙ (excitation voltage 175 V), (e) the fundamental wave A₁, and (f) the high order harmonics Aₙ (excitation voltage 200 V) during the whole fatigue process.</figcaption></figure>

Generally, the experiment results can also be divided into three stages: (I) the early dislocation fatigue stage, (II) the micro-crack sprouting stage, and (III) the crack expansion stage, which is consistent with the theoretical model results, as shown in the above figures. While the loading cycle was about less than 10,000, harmonic components increased obviously with continuous loading, which can be treated as the early stage of fatigue. In this stage, acoustic nonlinearity mainly came from dislocations, and it can be used to assess the early stage fatigue degree, as explained in the theoretical model of Region I. Further loading led to micro-cracks sprouting inside of the metals, and the occurrence of a weak interface reduced the transmitted wave energy. Therefore, the increasing
trend of the harmonics amplitude became slow even with higher dislocation density, which was the transition stage. The boundary line of this stage was between 10,000 and 15,000 cycles, approximately, which corresponds to the middle part in Figure 4. While the loading cycle was about 20,000, cracks were seen clearly at the end of the preset grooving, which meant the beginning of the crack expansion stage. The amplitude of each order harmonic declined quickly during this stage, the results of which are consistent with the theoretical results of Region III in Figure 4.

## 4. Conclusions

Acoustic nonlinearity during the fatigue process in metals can be divided into three stages: (I) the early dislocation fatigue stage, (II) the micro-crack sprouting stage, and (III) the crack expansion stage, which can be described by nonlinear spring models combined with dislocation theory. Varying spring stiffness is chosen to connect different stages in the fatigue process. For the early dislocation fatigue stage, spring stiffness approaches infinity, and acoustic nonlinearity mainly results from dislocations, the results of which can be used for early fatigue evaluation. Once micro-cracks sprout in Region II, spring stiffness will drop, which reduces the transmitted wave energy and makes the harmonics fall. The decreasing trends become steep quickly, owing to the crack's expansion, and additional spectrum components are generated in the results with crack occurrence in Region III. These may be the effects of the scattering wave on rough crack interface or the spatial variation in the stiffness of the springs. This phenomenon is complicated, and studying its nonlinear acoustic mechanism is our next work.

**Author Contributions:** Conceptualization, W.L.; validation, X.W. and W.X.; experiments' design and sample preparation, W.L. and X.W.; data curation, W.L. and W.X.; original draft preparation, W.L.; review and editing, X.W. and W.X.; supervision, X.W. and W.X.; funding acquisition, X.W.

**Funding:** This research was funded by the National Natural Science Foundation of China (11674351) and China Scholarship Council (201704910691).

**Acknowledgments:** The authors would like to express their gratitude for the help with the fatigue tests by the Technical Institute of Physics and Chemistry, Chinese Academy of Sciences.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Klepka, A.; Staszewski, W.; Jenal, R.; Szwedo, M.; Iwaniec, J.; Uhl, T. Nonlinear acoustics for fatigue crack detection—Experimental investigations of vibro-acoustic wave modulations. *Struct. Health Monit.* **2012**, *11*, 197–211. [CrossRef]

2. Tsyfansky, S.; Beresnevich, V. Non-linear vibration method for detection of fatigue cracks in aircraft wings. *J. Sound Vib.* **2000**, *236*, 49–60. [CrossRef]

3. Palitsagar, S.; Das, S.; Parida, N.; Bhattacharya, D. Non-linear ultrasonic technique to assess fatigue damage in structural steel. *Scr. Mater.* **2006**, *55*, 199–202. [CrossRef]

4. Nagy, P.B. Fatigue damage assessment by nonlinear ultrasonic materials characterization. *Ultrasonics* **1998**, *36*, 375–381. [CrossRef]

5. Yang, Z.; Tian, Y.; Li, W.; Zhou, H.; Zhang, W. Experimental Investigation of the Acoustic Nonlinear Behavior in Granular Polymer Bonded Explosives with Progressive Fatigue Damage. *Materials* **2017**, *10*, 660. [CrossRef] [PubMed]

6. Landau, L.D.; Lifshitz, E.M. *Theory of Elasticity*, 3rd ed.; Pergamon Press: New York, NY, USA, 1959; pp. 104–107.

7. Murnaghan, F.D. *Finite Deformations of an Elastic Solid*; John Wiley: New York, NY, USA, 1951; pp. 60–66.

8. Granato, A.V.; Lücke, K. The Vibrating String Model of Dislocation Damping. In *Physical Acoustics*; Mason, W.P., Ed.; Academic Press: New York, NY, USA, 1966; Volume 3, pp. 225–276.

9. Cantrell, J.H.; Yost, W.T. Nonlinear ultrasonic characterization of fatigue microstructures. *Int. J. Fatigue* **2001**, *23*, 487–490. [CrossRef]

10. Cantrell, J.H.; Yost, W.T. Acoustic harmonic generation from fatigue-induced dislocation dipoles. *Philos. Mag. A* **1994**, *69*, 315–326. [CrossRef]
11. Mayergoyz, I.D. Hysteresis models from the mathematical and control theory points of view. *J. Appl. Phys.* **1985**, *57*, 3803–3805. [CrossRef]

12. Gonzalez-Valadez, M.; Baltazar, A.; Dwyer-Joyce, R. Study of interfacial stiffness ratio of a rough surface in contact using a spring model. *Wear* **2010**, *268*, 373–379. [CrossRef]

13. Wu, X.; Li, H.; Lv, W.; Zhang, S.; Zhang, D. Numerical analysis on acoustic nonlinear phenomena between rough contact interfaces. In Proceedings of the 2016 IEEE International Ultrasonics Symposium (IUS), Tours, France, 18–21 September 2016; pp. 1–4.

14. An, Z.; Wang, X.; Deng, M.; Mao, J.; Li, M. A nonlinear spring model for an interface between two solids. *Wave Motion* **2013**, *50*, 295–309. [CrossRef]

15. Takahashi, S.; Motegi, R. Measurement of third-order elastic constants and applications to loaded structural materials. *SpringerPlus* **2015**, *4*, 325. [CrossRef] [PubMed]

16. Richardson, J.M. Harmonic generation at an unbonded interface—I. Planar interface between semi-infinite elastic media. *Int. J. Eng. Sci.* **1979**, *17*, 73–85. [CrossRef]

17. Buck, O. Acoustic harmonic generation at unbonded interfaces and fatigue cracks. *J. Acoust. Soc. Am.* **1978**, *33*, 371–373.

18. Cantrell, J.H.; Pan, E. Substructural organization, dislocation plasticity and harmonic generation in cyclically stressed wavy slip metals. *Proc. R. Soc. Lond. A Math. Phys. Eng. Sci.* **2004**, *460*, 757–780. [CrossRef]

19. Cantrell, J.H.; Yost, W.T. Effect of precipitate coherency strains on acoustic harmonic generation. *J. Appl. Phys.* **1997**, *81*, 2957–2962. [CrossRef]

20. Deng, M. Cumulative second-harmonic generation of Lamb-mode propagation in a solid plate. *J. Appl. Phys.* **1999**, *85*, 3051–3058. [CrossRef]

21. Pei, J.; Deng, M. Assessment of accumulated fatigue damage in solid plates using nonlinear Lamb wave approach. *Appl. Phys. Lett.* **2007**, *90*, 121902.

22. Shah, A.; Ribakov, Y. Non-linear ultrasonic evaluation of damaged concrete based on higher order harmonic generation. *Mater. Des.* **2009**, *30*, 4095–4102. [CrossRef]

23. Yan-Xun, X.; Fu-Zhen, X.; Ming-Xi, D. Evaluation of Thermal Degradation Induced Material Damage Using Nonlinear Lamb Waves. *Chin. Phys. Lett.* **2010**, *27*, 016202. [CrossRef]

24. Khelladi, H.; Rahmi, F. Perturbation methods for the spectral analysis of a weakly nonlinear acoustic field generated by a transient insonation. *J. Acoust. Soc. Am.* **2013**, *19*, 045079.

25. Liu, G.; Jayathilake, P.G.; Khoo, B.C. Perturbation method for the second-order nonlinear effect of focused acoustic field around a scatterer in an ideal fluid. *Ultrasonics* **2014**, *54*, 576–585. [CrossRef] [PubMed]

26. Momani, S.; Odibat, Z. Homotopy perturbation method for nonlinear partial differential equations of fractional order. *Phys. Lett. A* **2007**, *365*, 345–350. [CrossRef]

27. Blackstock, D.T. Convergence of the Keck-Beyer Perturbation Solution for Plane Waves of Finite Amplitude in a Viscous Fluid. *J. Acoust. Soc. Am.* **1966**, *39*, 411–413. [CrossRef]

28. Born, M. Dynamical Theory of Crystal Lattices. *Am. J. Phys.* **1955**, *23*, 474. [CrossRef]

<figure><img src="image_15.png" /></figure>