Interaction of Lamb Wave Modes with Weak Material Nonlinearity: Generation of Symmetric Zero-Frequency Mode

Abstract

The symmetric zero-frequency mode induced by weak material nonlinearity during Lamb wave propagation is explored for the first time. We theoretically confirm that, unlike the second harmonic, phase-velocity matching is not required to generate the zero-frequency mode and its signal is stronger than those of the nonlinear harmonics conventionally used, for example, the second harmonic. Experimental and numerical verifications of this theoretical analysis are conducted for the primary S₀ mode wave propagating in an aluminum plate. The existence of a symmetric zero-frequency mode is of great significance, probably triggering a revolutionary progress in the field of non-destructive evaluation and structural health monitoring of the early-stage material nonlinearity based on the ultrasonic Lamb waves.

1. Introduction

Material non-destructive evaluation and structural health monitoring during the early stage of material degradation are crucial for structural integrity and safety [1]. The appearance of nonlinear effects in elastic wave propagation is one of the most notable and sensitive indicators for the onset of material damage [2,3,4,5,6]. Previous works have shown that the techniques based on nonlinear Lamb waves can be utilized for long-range and in-depth inspection [7,8,9] through characterizing acoustic nonlinearity, which in turn can be regarded as an indicator of the early-stage material nonlinearity caused by fatigue plasticity, precipitation hardening of thermal aging, corrosion, etc. [10,11,12,13,14]. In general, a nonlinear response is generated in the region where the early-stage material nonlinearity occurs and then several nonlinear wave components may appear including subharmonics or higher-order harmonics. All of these possible nonlinear components can be used to evaluate or monitor the material nonlinearity by exploiting their higher sensitivity to the changes of the material microstructure [15] compared with the conventional techniques based on linear waves [7,16]. In the research field of the nonlinear Lamb waves, the observation and utilization of the second harmonic is the focus [1,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]. However, an accurate experimental implementation of nonlinear Lamb waves is difficult due to the inherent multi-mode and dispersive nature. Also, the amplitude of the second harmonic of Lamb waves is cumulative only when certain conditions (phase-velocity matching, non-zero power flux or even group-velocity matching) are satisfied [7,25,27,28]. Furthermore, the intensity of these nonlinear wave components is generally low. Even when the signals are not submerged by noises, they are still difficult to be detected and estimated [2]. In fact, the propagation of an ultrasonic wave through a nonlinear material not only leads to the generation of the harmonics of the original wave but also the generation of a static displacement component (also called direct current component or zero-frequency mode) [29,30,31,32,33,34,35,36,37]. Narasimha et al. [34,35], Jacob et al. [36] and Nagy et al. [37] have shown that the zero-frequency mode varies linearly with the propagation distance. However, their studies are mainly confined to longitudinal acoustic waves.

To overcome these limitations, Sun et al. (theory) [38] and Wan et al. (simulation) [39] have shown that the signal of the zero-frequency mode is stronger than that of the traditional nonlinear harmonics for Lamb waves. Unlike the second harmonic, phase-velocity matching is not required for the zero-frequency mode accumulation. We theoretically, experimentally and numerically demonstrate the evidence, the efficiency and the feasibility of this concept in this paper.

The interaction of the elastic Lamb waves with a nonlinear plate medium is a fascinating and extremely complex process [5], during which a finite-amplitude wave propagates in an isotropic and homogeneous plate with stress-free surfaces. For simplicity, it is assumed that the wave propagates in the x₁-direction and the particle motions are in the x₁–x₃ plane (Figure 1a). Geometrical and weak material nonlinearities are theoretically considered here and the total displacement field is assumed to be the sum of a primary wave ${\mathbf{u}}_{(1)}$ (at frequency $\omega$) and a second harmonic wave ${\mathbf{u}}_{(2)}$ (at frequency 2 $\omega$) based on the perturbation condition $\Vert u_{(2)}\Vert \ll \Vert u_{(1)}\Vert$ [15,21,22,40,41]. The perturbation method is used to solve the nonlinear motion equations for the harmonic generation in the waveguide. Solutions for the second harmonic, sum- and difference-frequency components are obtained via modal decomposition [28]. This procedure originally assumes that the nonlinearly generated secondary wave fields (perturbation solution) can be expressed as a superposition of the Lamb wave modes, since an orthogonality exists between different Lamb wave modes (completeness of Lamb wave modes is also assumed). It is obvious that the zero-frequency mode is orthogonal to other Lamb wave modes. Only when the zero-frequency mode is considered, the completeness of the mode expansion for the Lamb wave modes can be achieved. However, previous studies [15,21,25,28,40,41] ignored the solution for the zero-frequency mode theoretically. The present work firstly investigates the derivation of sum-frequency, difference-frequency, second harmonic and zero-frequency components of nonlinear Lamb waves in a plate with weak material nonlinearity. Then, the corresponding experimental and simulation results are demonstrated.

2. Formulations

2.1. Nonlinear Wave Equation

Following Lima et al. [21,28,40,41], the Green-Lagrange strain tensor, $\mathbf{E}$, is related to the displacement tensor, $\mathbf{u}(\mathbf{X},t)$, by

$$E=\frac{1}{2}(H+H^T+H^{T}H),H=\nabla u$$

(1)

where, $\mathbf{X}$ denotes the position of the material points in the reference configuration and $\mathbf{H}$ is the displacement gradient. Note that a black and bold letter denotes a tensor in this paper. The following strain energy, $W(\mathbf{E})$, proposed by Landau and Lifshitz [42] for a nonlinear hyperelastic isotropic solid, is employed

$$W(E)=\frac{1}{2}\lambda {[tr(\mathbf{E})]}^2+\mu tr(E^2)+\frac{1}{3}\mathcal{C}{[tr(E)]}^3+\mathcal{B}tr(E)tr(E^2)+\frac{1}{3}\mathcal{A}tr(E^3)$$

(2)

where, $\lambda$ and $\mu$ are Lame’s constants, $\mathcal{A}$, $\mathcal{B}$ and $\mathcal{C}$ denote the third-order elastic constants, which is related to material nonlinearity and $tr[·]$ represents the trace of the bracketed tensor. The second Piola-Kirchhoff stress tensor, $\mathbf{S}$, can be derived in terms of Lagrangian strain tensor, E, as

$$S=\frac{\partial W(E)}{\partial E}=\lambda tr(E)I+2\mu E+\mathcal{C}{[tr(E)]}^{2}I+\mathcal{B}tr(E^2)I+2\mathcal{B}tr(E)E+\mathcal{A}{E}^2$$

(3)

where, $\mathbf{I}$ is the identity tensor. Substitution of Equation (1) into (3) yields the second-order form of $\mathbf{H}$,

$$\begin{array}{c}S=\frac{\lambda }{2}tr(H+H^T)I+\mu (H+H^T)+\frac{\lambda }{2}tr(H^{T}H)I+\mu H^{T}H+\mathcal{C}{[tr(H)]}^{2}I+\\ \frac{\mathcal{B}}{2}tr(H^2+H^{T}H)I+\mathcal{B}tr(H)(H+H^T)+\frac{\mathcal{A}}{4}(H^2+H^{T}H+H{H}^T+{(H^T)}^2)\end{array}$$

(4)

$\mathbf{S}$ can be decomposed into linear ($S^L$) and nonlinear ($S^{NL}$) components by

$$\begin{array}{c}S=S^L+S^{NL},S^L=\lambda tr(H)I+\mu (H+H^T),\hfill \\ S^{NL}=\frac{\lambda }{2}tr(H^{T}H)I+\mu H^{T}H+\mathcal{C}{[tr(H)]}^{2}I+\frac{\mathcal{B}}{2}tr(H^2+H^{T}H)I\hfill \\ +\mathcal{B}tr(H)(H+H^T)+\frac{\mathcal{A}}{4}(H^2+H^{T}H+H{H}^T+{(H^T)}^2).\hfill \end{array}$$

(5)

The equation of the wave motion is more conveniently written in terms of the first Piola-Kirchhoff stress tensor, $\mathbf{T}$, which is related to the second Piola-Kirchhoff stress tensor, $\mathbf{S}$, through the deformation gradient tensor: $\mathbf{F}=\mathbf{H}+\mathbf{I}$, by

$$T=F\cdot S=F\cdot (S^L+S^{NL})=(H+I)\cdot (S^L+S^{NL})$$

(6)

Then, $\mathbf{T}$ can also be decomposed into linear ($T^L(\mathbf{H})$) and nonlinear ($T^{NL}(\mathbf{H})$) components by

$$\begin{array}{c}T(H)=T^L(H)+T^{NL}(H),T^L(H)=S^L=\lambda tr(H)I+\mu (H+H^T),\hfill \\ T^{NL}(H)=H{S}^L+S^{NL}=\frac{\lambda }{2}tr(H^{T}H)I+\mathcal{C}{[tr(H)]}^{2}I+\mathcal{B}tr(H)H^T+\frac{\mathcal{A}}{4}{(H^T)}^2\hfill \\ +\frac{B}{2}tr(H^2+H^{T}H)I+(\lambda +\mathcal{B})tr(H)H+(\mu +\frac{\mathcal{A}}{4})(H^2+H^{T}H+H{H}^T)\hfill \end{array}$$

(7)

where, only the terms up to the second-order are kept.

Considering the wave motion in an isotropic homogeneous plate (Figure 1a), the wave motion equation and stress-free boundary condition can be written as

$$\nabla \cdot [T(H)]=\rho \ddot{u},{[T(H)\cdot n_3]|}_{x_3=\pm h}=0$$

(8)

where, $n_3$ is a unit vector parallel to the $x_3$ axis and $2h$ is the thickness of the plate as shown in Figure 1a.

Under the condition of weak nonlinearity, Equation (8) can be solved using a perturbation method through decomposing the wave field into the primary and secondary components

$$u(X,t)=u_{(1)}(X,t)+u_{(2)}(X,t),\Vert u_{(2)}(X,t)\Vert \ll \Vert u_{(1)}(X,t)\Vert$$

(9)

where, $u_{(1)}(\mathbf{X},t)$ and $u_{(2)}(\mathbf{X},t)$ represent the primary and secondary wave field, respectively and $\Vert ·\Vert$ represents the amplitude of the corresponding wave field.

When substituting Equations (1), (7) and (9) into Equation (8), we obtain

$$\begin{array}{c}(\lambda +2\mu )\nabla \nabla \cdot [u_{(1)}(X,t)+u_{(2)}(X,t)]-\mu \nabla \times \nabla \times [u_{(1)}(X,t)\hfill \\ +u_{(2)}(X,t)]+\nabla \cdot [T^{NL}(H)]=\rho [{\ddot{u}}_{(1)}(X,t)+{\ddot{u}}_{(2)}(X,t)],\hfill \\ {[T(H)\cdot n_3]|}_{x_3=\pm h}={[T^L(H)\cdot n_3+T^{NL}(H)\cdot n_3]|}_{x_3=\pm h}=0.\hfill \end{array}$$

(10)

Now, Equation (10) can be divided into two parts, that is, the first-order and second-order approximation parts. The first-order approximation is

$$\begin{array}{c}(\lambda +2\mu )\nabla \nabla \cdot u_{(1)}(X,t)-\mu \nabla \times \nabla \times u_{(1)}(X,t)=\rho {\ddot{u}}_{(1)}(X,t),\hfill \\ {T^L[H(u_{(1)})]\cdot n_3|}_{x_3=\pm h}=0\hfill \end{array}$$

(11)

which gives the primary wave field and the solution is

$$\begin{array}{cc}\hfill \mathrm{SH}\mathrm{mode}:& q=\frac{n\pi }{2h},\{\begin{array}{c}n\in \left\{0,2,4,\dots \right\}\mathrm{symmteric},\hfill \\ n\in \left\{1,3,5,\dots \right\}\mathrm{antisymmetric},\hfill \end{array}\hfill \\ \mathrm{RL}\mathrm{mode}:\hfill & \frac{\tan (qh)}{\tan (ph)}=-{[\frac{4k^{2}pq}{{(k^2-q^2)}^2}]}^{\pm 1},\{\begin{array}{c}+1,\mathrm{symmteric},\hfill \\ -1,\mathrm{antisymmteric},\hfill \end{array}\hfill \\ & p^2=\frac{{\omega }^2}{c_L^2}-k^2\mathrm{and}{q}^2=\frac{{\omega }^2}{c_T^2}-k^2\hfill \end{array}$$

(12)

where, $c_L$ and $c_T$ are the velocities of the elastic longitudinal and transverse waves, respectively, SH and RL represent shear-horizontal and Rayleigh-Lamb, respectively. The second-order approximation is

$$\begin{array}{c}(\lambda +2\mu )\nabla \nabla \cdot u_{(2)}(X,t)-\mu \nabla \times \nabla \times u_{(2)}(X,t)-\rho {\ddot{u}}_{(2)}(X,t)=-\nabla \cdot \left\{T^{NL}[H(u_{(1)},u_{(1)})]\right\},\\ {T^L[H(u_{(2)})]\cdot n_3=-T^{NL}[H(u_{(1)},u_{(1)})]\cdot n_3]|}_{x_3=\pm h}\end{array}$$

(13)

$T^{NL}[\mathbf{H}(u_{(1)},u_{(1)})]$ and $\nabla \cdot \left\{T^{NL}[H(u_{(1)},u_{(1)})]\right\}$ which contain quadratic terms in $u_{(1)}$ can be taken as the force terms. Note that the higher-order force terms $T^{NL}[\mathbf{H}(u_{(1)},u_{(2)})]$, $T^{NL}[\mathbf{H}(u_{(2)},u_{(2)})]$, $\nabla \cdot \left\{T^{NL}[H(u_{(1)},u_{(2)})]\right\}$ and $\nabla \cdot \left\{T^{NL}[H(u_{(2)},u_{(2)})]\right\}$ are ignored. Then, Equation (13) can be regarded as a forced waveguide equation and be written as

$$\begin{array}{c}(\lambda +2\mu )\nabla \nabla \cdot u_{(2)}(X,t)-\mu \nabla \times \nabla \times u_{(2)}(X,t)-\rho {\ddot{u}}_{(2)}(X,t)=-f^{(1,1)},\\ {{S^{(2)}\cdot n_3|}_{x_3=\pm h}=-S^{(1,1)}\cdot n_3|}_{x_3=\pm h}.\end{array}$$

(14)

where, $S^{(2)}=T^L[\mathbf{H}(u_{(2)})]$, $S^{(1,1)}=T^{NL}[\mathbf{H}(u_{(1)},u_{(1)})]$ and $f^{(1,1)}=\nabla \cdot \left\{T^{NL}[H(u_{(1)},u_{(1)})]\right\}$. The expressions for $S^{(1,1)}$ and $f^{(1,1)}$ can be derived by Equations (1) and (7) as

$$\begin{array}{l}{S}_{ii}^{(1,1)}=\frac{\lambda }{2}\frac{\partial u_{(1)k}}{\partial x_l}\frac{\partial u_{(1)k}}{\partial x_l}+\mathcal{C}\frac{\partial u_{(1)k}}{\partial x_k}\frac{\partial u_{(1)l}}{\partial x_l}+\mathcal{B}\frac{\partial u_{(1)k}}{\partial x_k}\frac{\partial u_{(1)i}}{\partial x_i}+\\ \frac{\mathcal{B}}{2}(\frac{\partial u_{(1)k}}{\partial x_l}\frac{\partial u_{(1)k}}{\partial x_l}+\frac{\partial u_{(1)k}}{\partial x_l}\frac{\partial u_{(1)l}}{\partial x_k})+\frac{\mathcal{A}}{4}\frac{\partial u_{(1)i}}{\partial x_k}\frac{\partial u_{(1)k}}{\partial x_i}+(\lambda +\mathcal{B})\frac{\partial u_{(1)k}}{\partial x_k}\frac{\partial u_{(1)i}}{\partial x_i}+\\ (\mu +\frac{\mathcal{A}}{4})(\frac{\partial u_{(1)i}}{\partial x_k}\frac{\partial u_{(1)i}}{\partial x_k}+\frac{\partial u_{(1)k}}{\partial x_i}\frac{\partial u_{(1)k}}{\partial x_i}+\frac{\partial u_{(1)i}}{\partial x_k}\frac{\partial u_{(1)k}}{\partial x_i})+o(u_{(1)i}{}^3),\\ S_{ij}^{(1,1)}=\mathcal{B}\frac{\partial u_{(1)k}}{\partial x_k}\frac{\partial u_{(1)j}}{\partial x_i}+\frac{\mathcal{A}}{4}\frac{\partial u_{(1)j}}{\partial x_k}\frac{\partial u_{(1)k}}{\partial x_i}+(\lambda +\mathcal{B})\frac{\partial u_{(1)k}}{\partial x_k}\frac{\partial u_{(1)i}}{\partial x_j}+\\ (\mu +\frac{\mathcal{A}}{4})(\frac{\partial u_{(1)i}}{\partial x_k}\frac{\partial u_{(1)j}}{\partial x_k}+\frac{\partial u_{(1)k}}{\partial x_i}\frac{\partial u_{(1)k}}{\partial x_j}+\frac{\partial u_{(1)i}}{\partial x_k}\frac{\partial u_{(1)k}}{\partial x_j})+o(u_{(1)i}{}^3),\mathrm{for}i\ne j\end{array}$$

(15)

$$\begin{array}{l}{f}_i^{(1,1)}=(\mu +\frac{\mathcal{A}}{4})(\frac{{\partial }^2{u}_{(1)l}}{\partial x_k^2}\frac{\partial u_{(1)l}}{\partial x_i}+\frac{{\partial }^2{u}_{(1)l}}{\partial x_k^2}\frac{\partial u_{(1)i}}{\partial x_l}+2\frac{{\partial }^2{u}_{(1)i}}{\partial x_l\partial x_k}\frac{\partial u_{(1)k}}{\partial x_l})\\ +(\lambda +\mu +\frac{\mathcal{A}}{4}+\mathcal{B})(\frac{{\partial }^2{u}_{(1)l}}{\partial x_k\partial x_i}\frac{\partial u_{(1)l}}{\partial x_k}+\frac{{\partial }^2{u}_{(1)l}}{\partial x_l\partial x_k}\frac{\partial u_{(1)i}}{\partial x_k})+(\lambda +\mathcal{B})\frac{{\partial }^2{u}_{(1)i}}{\partial x_k^2}\frac{\partial u_{(1)l}}{\partial x_l}\\ +(\frac{\mathcal{A}}{4}+\mathcal{B})(\frac{{\partial }^2{u}_{(1)k}}{\partial x_l\partial x_k}\frac{\partial u_{(1)l}}{\partial x_i}+\frac{{\partial }^2{u}_{(1)l}}{\partial x_i\partial x_k}\frac{\partial u_{(1)k}}{\partial x_l})+(\mathcal{B}+2\mathcal{C})\frac{{\partial }^2{u}_{(1)l}}{\partial x_l\partial x_i}\frac{\partial u_{(1)k}}{\partial x_k}+o(u_{(1)i}^3)\end{array}$$

(16)

2.2. Second-Order Solution

2.2.1. Generation of Sum-Frequency, Difference-Frequency and Second Harmonic

Following de Lima et al. [21,28,40,41], consider two plane wave propagation modes with frequencies ${\omega }_a$ and ${\omega }_b$ and the corresponding wavenumbers $k_a$ and $k_b$, that are excited at $x_1=0$. The primary solution can be written as the form

$$u_{(1)}(x_1,x_3,t)=\frac{1}{2}{u}_{(a)}(x_3)e^{-i(k_a{x}_1-{\omega }_{a}t)}+\frac{1}{2}{u}_{(b)}(x_3)e^{-i(k_b{x}_1-{\omega }_{b}t)}+c.c.$$

(17)

where, $c.c.$ denotes the complex conjugate (Here, Equation (17) can also be written as $u_{(1)}=u_{(a)}(x_3)\cos (k_a{x}_1-{\omega }_{a}t)+u_{(b)}(x_3)\cos (k_b{x}_1-{\omega }_{b}t)$). Equation (17) can be introduced into Equations (15) and (16) to obtain

$$\begin{array}{l}{S}^{(1,1)}=S_{a+a}^{(1,1)}(x_3)e^{-i2(k_a{x}_1-{\omega }_{a}t)}+S_{b+b}^{(1,1)}(x_3)e^{-i2(k_b{x}_1-{\omega }_{b}t)}+\\ S_{a+b}^{(1,1)}(x_3)e^{-i[(k_a+k_b)x_1-({\omega }_a+{\omega }_b)t]}+S_{a-b}^{(1,1)}(x_3)e^{-i[(k_a-k_b)x_1-({\omega }_a-{\omega }_b)t]}+c.c.\\ f=f_{a+a}^{(1,1)}(x_3)e^{-i2(k_a{x}_1-{\omega }_{a}t)}+f_{b+b}^{(1,1)}(x_3)e^{-i2(k_b{x}_1-{\omega }_{b}t)}+\\ f_{a+b}^{(1,1)}(x_3)e^{-i[(k_a+k_b)x_1-({\omega }_a+{\omega }_b)t]}+f_{a-b}^{(1,1)}(x_3)e^{-i[(k_a-k_b)x_1-({\omega }_a-{\omega }_b)t]}+c.c.\end{array}$$

(18)

where, $S_{a+a}^{(1,1)}(x_3)$ and $f_{a+a}^{(1,1)}(x_3)$ are determined by the self-interaction of the excited mode $u_{(a)}(x_3)$, $S_{b+b}^{(1,1)}(x_3)$ and $f_{b+b}^{(1,1)}(x_3)$ are determined by the self-interaction of the excited mode $u_{(b)}(x_3)$ and $S_{a+b}^{(1,1)}(x_3)$, $S_{a-b}^{(1,1)}(x_3)$, $f_{a+b}^{(1,1)}(x_3)$ and $f_{a-b}^{(1,1)}(x_3)$ are due to the mutual interaction of $u_{(a)}(x_3)$ with $u_{(b)}(x_3)$. The terms with respect to $e^{-i2({\omega }_a+{\omega }_b)t}$ and $e^{-i2({\omega }_a-{\omega }_b)t}$ are the sum- and difference-frequency components, respectively and the terms related to $e^{-i2{\omega }_{a}t}$ and $e^{-i2{\omega }_{b}t}$ are the second harmonic components.

According to [28], to simplify the notation, Equation (18) can be written as the following form

$$\begin{array}{l}{S}^{(1,1)}=S_{\pm }^{(1,1)}(x_3)e^{-i[(k_a\pm k_b)x_1-({\omega }_a\pm {\omega }_b)t]}+\mathrm{c}.\mathrm{c}.,\\ f=f_{\pm }^{(1,1)}(x_3)e^{-i[(k_a\pm k_b)x_1-({\omega }_a\pm {\omega }_b)t]}+\mathrm{c}.\mathrm{c}.\end{array}$$

(19)

Following Auld [43] and de Lima and Hamilton [28], the secondary solution is derived as follows

$$\begin{array}{cc}\hfill v_{(2)}(X,t)& \hfill =\frac{\partial u_{(2)}(X,t)}{\partial t}=\frac{\partial }{\partial t}[\frac{1}{2}{\displaystyle \sum _{m=1}^{\infty }{A}_m}(x_1)u_m(x_3)e^{-2i({\omega }_a\pm {\omega }_b)t}+c.c.],\\ & =\frac{1}{2}{\displaystyle \sum _{m=1}^{\infty }{a}_m}(x_1)v_m(x_3)e^{-2i({\omega }_a\pm {\omega }_b)t}+c.c.,\hfill \\ & 4P_{mn}[\frac{d}{d{x}_1}+i{k}_n^{*}]a_m(x_1)=f_n{e}^{-i(k_a\pm k_b)x_1}\hfill \end{array}$$

(20)

where,

$$\begin{array}{c}{a}_m(x_1)=\frac{f_n}{4P_{mn}}\{\begin{array}{cc}\frac{i}{[(k_a\pm k_b)-k_n^{*}]}[e^{-i(k_a\pm k_b)x_1}-e^{-i{k}_n^{*}{x}_1}],\hfill & \mathrm{for}{k}_n^{*}\ne k_a\pm k_b\hfill \\ x_1{e}^{-i(k_a\pm k_b)x_1},\hfill & \mathrm{for}{k}_n^{*}=k_a\pm k_b\hfill \end{array}\\ P_{mn}\equiv -\frac{1}{4}{\displaystyle {\int }_{-h}^{+h}(\frac{v_n^{*}}{2}\cdot \frac{T_m}{2}+\frac{v_m}{2}\cdot \frac{T_n^{*}}{2}\cdot )\cdot {\widehat{x}}_{1}d{\widehat{x}}_3},P_{mn}=0\mathrm{if}{k}_m\ne k_n^{*},\\ f_n\equiv f_n^{surfce}+f_n^{volume},f_n^{surfce}\equiv \frac{1}{2}{[v_n^{*}(x_3)\cdot S_{\pm }^{(1,1)}(x_3)]\cdot {\widehat{x}}_3|}_{x_3=-h}^{x_3=+h},\\ f_n^{volume}\equiv -\frac{1}{2}{\displaystyle {\int }_{-h}^{+h}{v}_n^{*}(x_3)\cdot f_{\pm }^{(1,1)}(x_3)}d{\widehat{x}}_3\\ T_m\equiv \lambda \mathrm{tr}[\nabla u_m(x_3)]I+\mu (\nabla u_m(x_3)+u_m(x_3)\nabla )\end{array}$$

(21)

The subscript $m$ and $n$ represent the $m^{\mathrm{th}}$ and $n^{\mathrm{th}}$ modes at the frequency ${\omega }_a\pm {\omega }_b$, the superscript $*$ denotes the complex conjugate of the corresponding physical variables, $v_m$ (short for $v_m(x_3)$) is the particle velocity of the $m^{\mathrm{th}}$ mode and $T_m$, $\mathbf{u}$, $S^{(1,1)}$ and $f^{1,1)}$ can also be expressed as a linear combination of the corresponding waveguide modes (not presented here) similar to that of $\mathbf{v}$.

According to the above equation $P_{mn}=0\mathrm{if}{k}_m\ne k_n^{*},$ a propagating mode $m$ is orthogonal to all modes except itself and $P_{mm}$ is the complex power flux of the $m^{\mathrm{th}}$ propagating mode in the direction ${\widehat{x}}_3$. $f_n^{surfce}$ and $f_n^{volume}$ are interpreted as the power flux through the surface and the volume, respectively, due to the primary wave. When $k_n^{*}=k_a\pm k_b$ and $f_n\ne 0$, the amplitude of the secondary solution increases linearly in the direction of the propagation. Hence, two conditions must be satisfied for the internal resonance: (1) phase-velocity matching, $k_n^{*}=k_a\pm k_b$; (2) non-zero power flux, $f_n\ne 0$.

On the basis of Equation (20), the sum- and difference-frequency components can be generated through internal resonance and the second harmonic generation can be considered as a special case of the sum-frequency generation, in which only a single mode is excited. Then the solution for the second harmonic can be written in the normal mode expansion form [15,22,25].

$$\begin{array}{c}{v}_{(2)}(x_1,x_3,t)={\displaystyle {\sum }_{n=1}^{\infty }{a}_n(x_1)}{v}_n(x_3)e^{-i2\omega t}+c.c.\\ a_n(x_1)=\frac{(f_n^{surface}+f_n^{volume})}{(4P_{nn})}\{\begin{array}{ccc}i[e^{-i2k{x}_1}-e^{-i{k}_n^{*}{x}_1}]/(k_n^{*}-2k)\hfill & \mathrm{for}& k_n^{*}\ne 2k\hfill \\ x_1{e}^{-i2k{x}_1}\hfill & \mathrm{for}& k_n^{*}=2k\hfill \end{array}\hfill \end{array}$$

(22)

where, $a_n(x_1)$ is the amplitude of the n^th second harmonic mode $v_n(x_3)$, $k$ and $k_n^{*}$ are the wave numbers of the primary mode and the n^th second harmonic mode, respectively, $P_{nn}$ is the power carried by the nth second harmonic mode and $f_n^{volume}$ and $f_n^{surface}$ are power fluxes from the primary mode to the second harmonic mode through the volume and surface of the plate, respectively. The detailed expressions for $f_n^{volume}$, $f_n^{surface}$ and $P_{nn}$ are complex and extremely lengthy and will not be presented here.

2.2.2. Symmetry Properties and Zero-Frequency Generation

According to [15], the following notations are introduced: $\mathrm{S}$ represents a unspecified, generic element of the set of symmetric functions in $x_3$, while $\mathrm{A}$ represents a generic, unspecified element of the set of antisymmetric functions in $x_3$. It is obvious that a derivation of a displacement component function with respect to $x_3$ changes the type of symmetry in $x_3$, while a derivation with respect to $x_1$ does not change the type of symmetry in $x_3$. Also, the following rules can be obtained easily:

$$\begin{array}{c}\mathrm{S}\cdot \mathrm{S}=\mathrm{S},\mathrm{S}\cdot \mathrm{A}=\mathrm{A}\cdot \mathrm{S}=\mathrm{A},\mathrm{A}\cdot \mathrm{A}=\mathrm{S},\\ \mathrm{S}+\mathrm{S}=\mathrm{S},\mathrm{A}+\mathrm{A}=\mathrm{A}\end{array}$$

(23)

Moreover, the parity of a symmetric mode can be represented as

$${\left\{\mathrm{S}0\mathrm{A}\right\}}^T$$

(24)

Now, all possible secondary modes can be considered by evaluating $f_n^{surface}$ and $f_n^{surface}$ based on parity. The parity of $\mathbf{H}$ can be written as

$$\begin{array}{c}\mathrm{SH}-\mathrm{S}\leftrightarrow \mathrm{symmetric}\mathrm{shear}-\mathrm{horizontal}\mathrm{mode},\mathrm{SH}-\mathrm{A}\leftrightarrow \mathrm{antisymmetric}\mathrm{shear}-\mathrm{horizontal}\mathrm{mode},\\ \mathrm{RL}-\mathrm{S}\leftrightarrow \mathrm{symmetric}\mathrm{Rayleigh}-\mathrm{Lamb}\mathrm{mode},\mathrm{RL}-\mathrm{A}\leftrightarrow \mathrm{antisymmetric}\mathrm{Rayleigh}-\mathrm{Lamb}\mathrm{mode},\\ H=H(\mathrm{SH}-\mathrm{S})+H(\mathrm{SH}-\mathrm{A})+H(\mathrm{RL}-\mathrm{S})+H(\mathrm{RL}-\mathrm{A}),\\ H(\mathrm{SH}-\mathrm{S})=\left[\begin{array}{ccc}\begin{array}{c}0\\ {\mathrm{S}}_a+{\mathrm{S}}_b+{\mathrm{S}}_a^{*}+{\mathrm{S}}_b^{*}\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}0\\ {\mathrm{A}}_a+{\mathrm{A}}_b+{\mathrm{A}}_a^{*}+{\mathrm{A}}_b^{*}\\ 0\end{array}\end{array}\right],H(\mathrm{SH}-\mathrm{A})=\left[\begin{array}{ccc}\begin{array}{c}0\\ {\mathrm{A}}_a+{\mathrm{A}}_b+{\mathrm{A}}_a^{*}+{\mathrm{A}}_b^{*}\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}0\\ {\mathrm{S}}_a+{\mathrm{S}}_b+{\mathrm{S}}_a^{*}+{\mathrm{S}}_b^{*}\\ 0\end{array}\end{array}\right],\\ H(\mathrm{RL}-\mathrm{S})=\left[\begin{array}{ccc}\begin{array}{c}{\mathrm{S}}_a+{\mathrm{S}}_b+{\mathrm{S}}_a^{*}+{\mathrm{S}}_b^{*}\\ 0\\ {\mathrm{A}}_a+{\mathrm{A}}_b+{\mathrm{A}}_a^{*}+{\mathrm{A}}_b^{*}\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}{\mathrm{A}}_a+{\mathrm{A}}_b+{\mathrm{A}}_a^{*}+{\mathrm{A}}_b^{*}\\ 0\\ {\mathrm{S}}_a+{\mathrm{S}}_b+{\mathrm{S}}_a^{*}+{\mathrm{S}}_b^{*}\end{array}\end{array}\right],H(\mathrm{RL}-\mathrm{A})=\left[\begin{array}{ccc}\begin{array}{c}{\mathrm{A}}_a+{\mathrm{A}}_b+{\mathrm{A}}_a^{*}+{\mathrm{A}}_b^{*}\\ 0\\ {\mathrm{S}}_a+{\mathrm{S}}_b+{\mathrm{S}}_a^{*}+{\mathrm{S}}_b^{*}\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}{\mathrm{S}}_a+{\mathrm{S}}_b+{\mathrm{S}}_a^{*}+{\mathrm{S}}_b^{*}\\ 0\\ {\mathrm{A}}_a+{\mathrm{A}}_b+{\mathrm{A}}_a^{*}+{\mathrm{A}}_b^{*}\end{array}\end{array}\right].\end{array}$$

(25)

Note that the complex conjugate components are considered here. Hence the parity of $S^{(1,1)}(x_3)$ and $f^{(1,1)}(x_3)$ is determined as follows

(26)

Then, $f_n^{surface}$ and $f_n^{surface}$ can be evaluated and only two terms below holding non-zero have been obtained

$$f_n^{surface}(\mathrm{RL}-\mathrm{S})={\frac{1}{2}\left\{\begin{array}{c}\mathrm{S}({\mathrm{S}}_a+{\mathrm{S}}_b+{\mathrm{S}}_a^{*}+{\mathrm{S}}_b^{*})\cdot ({\mathrm{A}}_a+{\mathrm{A}}_b+{\mathrm{A}}_a^{*}+{\mathrm{A}}_b^{*})+\mathrm{A}\cdot \\ [{({\mathrm{A}}_a+{\mathrm{A}}_b+{\mathrm{A}}_a^{*}+{\mathrm{A}}_b^{*})}^2+({\mathrm{S}}_a+{\mathrm{S}}_b+{\mathrm{S}}_a^{*}+{\mathrm{S}}_b^{*}{)}^2]\end{array}\right\}|}_{x_3=-h}^{x_3=+h}={\frac{1}{2}\mathrm{A}|}_{x_3=-h}^{x_3=+h}\ne 0,$$

(27a)

$$\begin{array}{c}{f}_n^{volume}(\mathrm{RL}-\mathrm{S})=-\frac{1}{2}{\displaystyle {\int }_{-h}^{+h}\left\{\begin{array}{c}\mathrm{S}\left\{\begin{array}{c}{({\mathrm{A}}_a+{\mathrm{A}}_b+{\mathrm{A}}_a^{*}+{\mathrm{A}}_b^{*})}^2+({\mathrm{S}}_a+{\mathrm{S}}_b+{\mathrm{S}}_a^{*}+{\mathrm{S}}_b^{*}{)}^2+\\ \frac{\partial }{x_3}[({\mathrm{S}}_a+{\mathrm{S}}_b+{\mathrm{S}}_a^{*}+{\mathrm{S}}_b^{*})\cdot ({\mathrm{A}}_a+{\mathrm{A}}_b+{\mathrm{A}}_a^{*}+{\mathrm{A}}_b^{*})]\end{array}\right\}+\\ \mathrm{A}\left\{\begin{array}{c}({\mathrm{S}}_a+{\mathrm{S}}_b+{\mathrm{S}}_a^{*}+{\mathrm{S}}_b^{*})\cdot ({\mathrm{A}}_a+{\mathrm{A}}_b+{\mathrm{A}}_a^{*}+A_b^{*})+\frac{\partial }{\partial x_3}\\ [{({\mathrm{A}}_a+{\mathrm{A}}_b+{\mathrm{A}}_a^{*}+{\mathrm{A}}_b^{*})}^2+({\mathrm{S}}_a+{\mathrm{S}}_b+{\mathrm{S}}_a^{*}+{\mathrm{S}}_b^{*}{)}^2]\end{array}\right\}\end{array}\right\}}d{\widehat{x}}_3\\ =-\frac{1}{2}{\displaystyle {\int }_{-h}^{+h}\mathrm{S}}d{x}_3\ne 0\end{array}$$

(27b)

Summarizing Equation (27), it is concluded that both a symmetric and an antisymmetric primary mode can excite a symmetric secondary mode and a symmetric zero-frequency mode. In contrast, neither a symmetric nor an antisymmetric primary mode can excite an antisymmetric secondary mode and an antisymmetric zero-frequency mode. This conclusion is similar to that of a previous study [15]. However, their conclusions are limited to the second harmonic generation only. Here, from Equation (27), we can see that ${({\mathrm{S}}_a)}^2$, ${({\mathrm{S}}_b)}^2$, ${({\mathrm{S}}_a^{*})}^2$, ${({\mathrm{S}}_b^{*})}^2$, ${({\mathrm{A}}_a)}^2$, ${({\mathrm{A}}_b)}^2$, ${({\mathrm{A}}_a^{*})}^2$, ${({\mathrm{A}}_b^{*})}^2$, ${\mathrm{S}}_a{\mathrm{A}}_a$, ${\mathrm{S}}_b{\mathrm{A}}_b$, ${\mathrm{S}}_a^{*}{\mathrm{A}}_a^{*}$ and ${\mathrm{S}}_b^{*}{\mathrm{A}}_b^{*}$ drive the second harmonic frequency; ${\mathrm{S}}_a{\mathrm{S}}_b$, ${\mathrm{S}}_a{\mathrm{A}}_b$, ${\mathrm{S}}_b{\mathrm{A}}_a$, ${\mathrm{A}}_a{\mathrm{A}}_b$, ${\mathrm{S}}_a^{*}{\mathrm{S}}_b^{*}$, ${\mathrm{S}}_a^{*}{\mathrm{A}}_b^{*}$, ${\mathrm{A}}_a^{*}{\mathrm{S}}_b^{*}$ and ${\mathrm{A}}_a^{*}{\mathrm{A}}_b^{*}$ drive the sum-frequency; ${\mathrm{S}}_a{\mathrm{S}}_b^{*}$, ${\mathrm{S}}_a{\mathrm{A}}_b^{*}$, ${\mathrm{S}}_b{\mathrm{A}}_a^{*}$, ${\mathrm{S}}_b{\mathrm{S}}_a^{*}$, ${\mathrm{A}}_a{\mathrm{A}}_b^{*}$, ${\mathrm{A}}_b{\mathrm{A}}_a^{*}$, ${\mathrm{A}}_b{\mathrm{A}}_a^{*}$ and ${\mathrm{A}}_b{\mathrm{S}}_a^{*}$ drive the difference-frequency; ${\mathrm{S}}_a{\mathrm{S}}_a^{*}$, ${\mathrm{S}}_a{\mathrm{A}}_a^{*}$, ${\mathrm{S}}_b{\mathrm{S}}_b^{*}$, ${\mathrm{S}}_b{\mathrm{A}}_b^{*}$, ${\mathrm{A}}_a{\mathrm{A}}_a^{*}$, ${\mathrm{A}}_a{\mathrm{S}}_a^{*}$, ${\mathrm{A}}_b{\mathrm{A}}_b^{*}$ and ${\mathrm{S}}_b^{*}{\mathrm{A}}_b$ drive the zero-frequency.

2.2.3. Zero-Frequency Generation for a Primary Wave of a Single Mode

In order to obtain the solution of zero-frequency mode, one single mode is considered for simplicity. We write the primary displacement wave filed at frequency $\omega$:

$$u_{(1)}(x_1,x_3,t)=\frac{1}{2}{u}_{(1)}(x_3)e^{-i(k{x}_1-\omega t)}+\frac{1}{2}{u}_{(1)}(x_3)e^{i(k{x}_1-\omega t)}$$

(28)

Here, Equation (28) can also be written as $u_{(1)}(x_1,x_3,t)=u_{(1)}(x_3)\cos (k{x}_1-\omega t)$.

Hence, $\mathbf{H}$ has the form

$$H=\frac{1}{2}\left[\begin{array}{ccc}\begin{array}{c}-ik{u}_{(1)}{}_1(x_3)\\ -ik{u}_{(1)2}(x_3)\\ -ik{u}_{(1)3}(x_3)\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}\frac{d{u}_{(1)1}(x_3)}{d{x}_3}\\ \frac{d{u}_{(1)2}(x_3)}{d{x}_3}\\ \frac{d{u}_{(1)3}(x_3)}{d{x}_3}\end{array}\end{array}\right]e^{-i(k{x}_1-\omega t)}+\frac{1}{2}\left[\begin{array}{ccc}\begin{array}{c}ik{u}_{(1)1}(x_3)\\ ik{u}_{(1)2}(x_3)\\ ik{u}_{(1)3}(x_3)\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}\frac{d{u}_{(1)1}(x_3)}{d{x}_3}\\ \frac{d{u}_{(1)2}(x_3)}{d{x}_3}\\ \frac{d{u}_{(1)3}(x_3)}{d{x}_3}\end{array}\end{array}\right]e^{i(k{x}_1-\omega t)}$$

(29)

According to the procedure of Section 2.2.2, after ignoring the terms of the second harmonic components and setting $k_a=k_b=k,k_n^{*}=0,$ we obtain $P_{nn}^0$, $f_n^{0surface}$ and $f_n^{0volume}$ for Rayleigh-Lamb ($\mathrm{RL}$) mode wave as

$$P_{nn}^0\equiv -\frac{1}{4}{\displaystyle {\int }_{-h}^{+h}(\frac{v_n^{*}}{2}\cdot \frac{T_n}{2}+\frac{v_n}{2}\cdot \frac{T_n^{*}}{2}\cdot )\cdot {\widehat{x}}_{1}d{\widehat{x}}_3}=-\frac{1}{8}{\displaystyle {\int }_{-h}^{+h}({\mathrm{v}}_{\mathrm{n}i}{\mathrm{T}}_{ni})d{\widehat{x}}_3},$$

(30)

$$f_n^{0surface}={\left\{\begin{array}{l}\frac{1}{2}{v}_{n1}(x_3)[(\frac{2B+A+2\mu }{4})k^2{u}_{(1)1}{u}_{(1)3}+(\frac{4\lambda +4B+8\mu +A}{8})u_{(1)1,3}{u}_{(1)3,3}]+\\ v_{n3}(x_3)[(\frac{\lambda +2\mathcal{C}+2B}{4}){[k{u}_{(1)1}]}^2+(\frac{2\lambda +2B+4\mu +A}{8})\left\{{[k{u}_{(1)3}]}^2+{[u_{(1)1,3}]}^2\right\}\\ +(\frac{3\lambda +2\mathcal{C}+6\mu +6B+2A}{4}){[u_{(1)3,3}]}^2\end{array}\right\}|}_{x_3=-h}^{x_3=+h},$$

(31)

$$f_n^{0volume}=-\frac{1}{2}{\displaystyle {\int }_{-h}^{+h}\left\{\begin{array}{l}{v}_{n1}(x_3)\left\{\begin{array}{l}(\frac{2B+2\mu +A}{4})k^2[u_{(1)1,3}{u}_{(1)3}+u_{(1)1}{u}_{(1)3,3}]\\ (\frac{4\lambda +4B+8\mu +A}{8})[u_{(1)1,33}{u}_{(1)3,3}+u_{(1)1,3}{u}_{(1)3,33}]\end{array}\right\}+\\ v_{n3}(x_3)\left\{\begin{array}{l}(\frac{\lambda +2\mathcal{C}+2B}{2})[k{u}_{(1)1}{u}_{(1)1,3}]+(\frac{2\lambda +2B+4\mu +A}{4})\\ \left\{\begin{array}{c}{k}^2{u}_{(1)3}{u}_{(1)3,3}\\ +u_{(1)1,3}{u}_{(1)1,33}\end{array}\right\}+(\frac{3\lambda +2\mathcal{C}+6\mu +6B+2A}{2})[u_{(1)3,33}]\end{array}\right\}\end{array}\right\}d}{\widehat{x}}_3$$

(32)

where, for $\mathrm{RL}-\mathrm{S}$ mode,

$$\begin{array}{l}{v}_1=-i\omega [ik{A}_{2}cos(p{x}_3)+B_{1}q\cos (q{x}_3)]\\ v_3=-i\omega [-A_{2}p\sin (p{x}_3)-ik{B}_{1}sin(q{x}_3)]\\ {\mathrm{T}}_1=\mu [-2ik{A}_{2}p\sin (p{x}_3)+(k^2-q^2)B_{1}sin(q{x}_3)]\\ {\mathrm{T}}_3=[-\lambda (k^2+p^2)A_{2}cos(p{x}_3)-2\mu (p^2{A}_{2}cos(p{x}_3)+ikq{B}_1\cos (q{x}_3))]\end{array}\}$$

(33)

$$\begin{array}{l}{u}_{(1)1}=\left\{[ik{A}_{2}cos(p{x}_3)+B_{1}q\cos (q{x}_3)]\right\}\\ u_{(1)3}=\left\{[-A_{2}p\sin (p{x}_3)-ik{B}_{1}sin(q{x}_3)]\right\}\end{array}\}$$

(34)

and for $\mathrm{RL}-\mathrm{A}$ mode,

$$\begin{array}{l}{v}_1=-i\omega [ik{A}_{1}sin(p{x}_3)-B_{2}q\sin (q{x}_3)]\\ v_3=-i\omega [A_{1}p\cos (p{x}_3)-ik{B}_{2}cos(q{x}_3)]\\ {\mathrm{T}}_1=\mu [2ik{A}_{1}p\cos (p{x}_3)+(k^2-q^2)B_{2}cos(q{x}_3)]\\ {\mathrm{T}}_3=[-\lambda (k^2+p^2)A_{1}sin(p{x}_3)-2\mu (p^2{A}_{1}sin(p{x}_3)-ikq{B}_2\sin (q{x}_3))]\end{array}\}$$

(35)

$$\begin{array}{l}{u}_{(1)1}=\left\{[ik{A}_{1}sin(p{x}_3)-B_{2}q\sin (q{x}_3)]\right\}\\ u_{(1)3}=\left\{[A_{1}p\cos (p{x}_3)-ik{B}_{2}cos(q{x}_3)]\right\}\end{array}\}$$

(36)

where, $A_1$, $A_2$, $B_1$ and $B_2$ are constants [44]. Hence, the amplitude of the zero-frequency mode is obtained as (by setting $k_a=k_b=k,k_n^{*}=0$ in Equation (20))

$$a_n^0(x_1)=(f_n^{0surface}+f_n^{0volume})x_1/4P_{nn}^0$$

(37)

where, $f_n^{0volume}$ and $f_n^{0surface}$ are power fluxes from the primary mode to the zero-frequency mode through the volume and surface of the plate, respectively and $P_{nn}^0$ is the power carried by the zero-frequency mode. However, the expressions for $P_{nn}^0$, $f_n^{0surface}$ and $f_n^{0volume}$ are complex and extremely lengthy and will not be presented here.

2.2.4. Zero-Frequency Mode versus the Second Harmonic

According to Equations (22) and (37), the amplitude of the zero-frequency mode can grow linearly in the direction of the wave propagation and the condition of phase-velocity matching ($k_n^{*}=2k$) is no longer required (or the condition of phase-velocity matching is always satisfied). Only one condition for the accumulation of the symmetric zero-frequency mode is quantified, that is, the non-zero power flux.

The corresponding explicit solutions of Equations (22) and (37) are very complex. Therefore, the relationship between the amplitudes of the zero-frequency mode and the second harmonic is implicitly evaluated here. From Equations (22) and (37), we can write a sum amplitude of the zero-frequency mode and the second harmonic as $A_n(x_1)=a_n^0(x_1)+a_n(x_1)$ $=$ $\mathrm{g}(H^{\mathrm{T}}H,H{H}^{\mathrm{T}})$. By considering a one-dimensional system in $x_1$ (simplified as $x$) and then, $H^{\mathrm{T}}H={(\frac{du}{dx})}^2$ and $H{H}^{\mathrm{T}}={(\frac{du}{dx})}^2$, so we have $A_n(x)=\mathrm{g}({(\frac{du}{dx})}^2)\propto {(\frac{du}{dx})}^2$. By only considering the second-order approximation, $\mathrm{g}$ is a linear function of ${(\frac{du}{dx})}^2$ by neglecting the high-order terms. An approximate estimation can be simply made by setting $u=\sin (kx+\omega t)$ and then ${(\frac{du}{dx})}^2=k^2{\cos }^2(kx+\omega t)=\frac{k^2}{2}+\frac{k^2}{2}\cos 2(kx+\omega t)$. The constant term, that is, $\frac{k^2}{2}$, is related to the zero-frequency mode and the other term, $\frac{k^2}{2}\cos 2(kx+\omega t)$, denotes the second harmonic. Although this coefficient of the two wave modes is equal, the energy carried by the zero-frequency mode is two times larger than that of the second harmonic via integral operation. If the higher harmonics generated in this process are considered, more energy will flow into the zero-frequency mode while less energy will flow to the second harmonic. Based on the above analysis, we can conclude that the signal intensity of the zero-frequency mode should be much stronger than that of the second harmonic.

3. Experiment

Figure 2 is a schematic of the experimental setup. The experimental study was carried out on a rectangular aluminum plate (2500 mm × 68 mm × 2 mm). Due to inevitable dislocation, micro-voids and multi-poles during the manufacture process of the aluminum plate, the specimen had weak material nonlinearity and no additional damage was introduced in the experiments. The Ritec SNAP system (RAM 5000) with a high power gated amplifier was used. A Hanning windowed tone burst signal of 200 kHz (the carrier frequency), 10 cycles and 100 Vp-p amplitude was input to the transmitting transducer with a center frequency of 250 kHz, which consisted of a Plexiglas wedge and a z-cut piezoelectric (LiNbO3) wafer. The carrier frequency 200 kHz is in the bandwidth of this transmitting transducer. The desired S0 Lamb wave can be induced by the refraction of the bulk longitudinal wave from the Plexiglas wedge into the aluminum plate. Note that, for simplicity, we focused on S0 mode working as the primary wave (Figure 1b). Guided waves travelling from the transmitter to the receiver (a Plexiglas wedge and a z-cut piezoelectric (LiNbO3) wafer) over a certain distance were measured. The center frequency of the receiver is 500 kHz. The received signal was amplified by both a preamplifier with −20 dB gain and a receiver amplifier with 32 dB gain.

Typical experimental excitation signal in time and frequency domain are shown in Figure 3a,b. Noise analysis will be shown in Section 5. And typical experimental signals in the time-domain at the propagation distance of 150, 350, 550 and 750 mm are shown in Figure 4b, respectively. The result of the short-time Fourier transform (STFT) corresponding to the wave packet at 150 mm within 0.1 ms is demonstrated in Figure 4c, from which the components of both zero-frequency and the fundamental frequency can be identified clearly, whereas the second harmonic is too weak to be observed clearly. The results of the fast Fourier transform (FFT) are shown in Figure 4d, corresponding to eight propagation distances. The second harmonic and the zero-frequency mode were generated. The second harmonic is too weak to be observed at all signal extraction points. Again, the intensity of the zero-frequency mode is much stronger than that of the second harmonic. As the wave propagation distance increases, the amplitude of the zero-frequency mode increases while the amplitude of the fundamental wave decreases. The cumulative effect of the second harmonic is not obvious due to no strict matching of phase-velocity (see Figure 1d). Then, the energy of the primary wave is mainly transferred to the zero-frequency mode during the wave propagation.

4. Simulation

Numerical simulations on the propagation of Lamb waves in an aluminum plate with 2 mm thickness were also carried out using a commercial FEM software, that is, Abaqus, without considering the damping effect. Similar to [45], the Landau and Lifshitz model of hyper-elasticity was adopted by a user subroutine VUMAT, as generally used for material definition. The material properties are shown in

Table 1. Material properties of Al.

ρ (kg/m3)	λ (GPa)	μ (GPa)	$\mathcal{A}\left(\mathbf{GPa}\right)$	$\mathcal{B}\left(\mathbf{GPa}\right)$	$\mathcal{C}\left(\mathbf{GPa}\right)$
2704	70.3	26.96	−4160	−1310	−1505

. The simulation model of the plate is set to 2400 mm long, which is long enough to eliminate the reflected wave. The fixed boundary condition is applied to the right end to eliminate effect of rigid displacement during the simulation and the upper and lower surfaces of the plate are free. Prescribed uniform displacement is actuated on the left end of the plate to excite the desired primary S0 mode Lamb waves. The actuating function of the excitation signal is described as: $x(t)=\frac{A}{2}\sin (2\pi ft)(1-\cos (2\pi \frac{f}{N}t)$, where, f (=200 kHz) is the central frequency, N (=10) is the number of sinusoidal cycles in a pulse and A (=0.0001 mm) is the amplitude of tone-burst. Nineteen signal detection points are considered with a distance of 150, 250, 350, 450, 550, 650, 750, 850, 950, 1050, 1150, 1250, 1350, 1450, 1550, 1650, 1750, 1850 and 1950 mm from the left of the model, respectively.

Figure 5a shows the simulation results, that is, the wave packets at the propagation distances of 150, 350, 550 and 750 mm, respectively. The insert of Figure 5a illustrates the total displacement field including all frequency components, indicating a typical symmetrical deformation mode of the plate during wave propagation. The STFT result of the wave packet at 150 mm within 0.1 ms is demonstrated in Figure 5b, which is consistent with the experimental data shown in Figure 4c. In Figure 5b, the components of the zero-frequency, fundamental frequency and second harmonic can be observed, while the intensity of the second harmonic is much lower than those of the zero-frequency and fundamental frequency. Figure 5c shows the results of STFT at wave propagation distance 250, 550 and 850 mm, from which the intensity growth of the zero-frequency and second harmonic components with the distance can be identified, indicating their intrinsic cumulative characteristic. The FFT results of the wave packets at eight propagation distances are shown in Figure 5d. The intensity growth of the zero-frequency and second harmonic components can still be observed. Again, we can see that the intensity of the zero-frequency mode is much higher than that of the second harmonic. The above numerical results are basically consistent with the experimental data shown in Figure 4, though the material and damping properties of the experimental plate are different from those in the numerical model. In fact, the stronger zero-frequency mode compared with the second harmonic was also obtained in [45] through a COMSOL simulation; nevertheless, it was ignored. The results of our simulation are consistent with those of a previous study [39]. Physically, the zero-frequency displacement mode represents the irreversible shift of a particle from its original equilibrium position. To qualitatively pick up the displacement components of the zero-frequency mode, the following operation was performed. For all particles in a section along the through-thickness direction (nodes in FEM simulation), when a complete wave packet passing through this section, we discretized the wave packet into many discrete points in the time-domain and made the summation of u₁ (along x₁-axis in Figure 1a) and u₂ (along x₃-axis in Figure 1a). For all sections along the x₁-axis, that is, the different wave propagation distances, we repeated the above operation and obtained the displacement fields for zero-frequency mode. The results of the displacement fields for u₁ and u₂ are shown in Figure 5e,f, indicating the symmetry of the zero-frequency mode. Moreover, u₁ and u₂ are accumulated as the propagation distance increases. This characteristic is schematically demonstrated in Figure 5e. We note in passing that, by varying the amplitude and the frequency of the excitation signal, the trends of the obtained results are basically the same, indicating the high reliability of the above observed phenomena.

5. Acoustic Nonlinearity Parameter

According to previous studies [7,8,9,10,12,13,14,18,19,20,23,24], the measured acoustic nonlinearity parameter of Lamb waves can be expressed as $\beta \propto A_2/A_1^2$, where, $A_1$ and $A_2$ are the measured amplitudes of the fundamental and the second harmonic signals, respectively. Here, the expressions ${\beta }_0\propto A_0/A_1^2$ and ${\beta }_2\propto A_2/A_1^2$ are adopted to measure the acoustic nonlinearity, where, $A_0$ is the measured amplitude of the zero-frequency mode. It can be seen from Figure 3b that the experimental excitation signal contains an obvious zero-frequency component and a weak second harmonic component. This is caused by electro-instrument noises such as circuit quadratic nonlinearity, the contact nonlinearity between transducers (transmitter and piezoelectric wafer sensors) and specimen. In order to reduce the influences of these noises, the following data processing method was adopted. We set the amplitudes of zero-frequency, fundamental frequency and second harmonic of excitation signal as $A_{e0}$, $A_{e1}$ and $A_{e2}$ respectively, as shown in Figure 3b. For each received signal, the amplitudes of zero-frequency, fundamental frequency and second harmonic were represented as $A_{r0}$, $A_{r1}$ and $A_{r2}$ respectively, as shown in Figure 3c. Then the amplitudes of zero-frequency (second harmonic) of the received signal were rescaled using the equation $A_{r0}^{*}=A_{r0}-A_{e0}{A}_{r1}/A_{e1}$ ($A_{r2}^{*}=A_{r2}-A_{e2}{A}_{r1}/A_{e1}$). Finally, the acoustic nonlinearity parameter could be calculated by the equation ${\beta }_0=A_{r0}^{*}/A_{r1}^2$ (${\beta }_2=A_{r2}^{*}/A_{r1}^2$). Note that the data processing method was not used for the simulation data. The experimental and numerical ${\beta }_0$ and ${\beta }_2$ are shown in Figure 6a,b, respectively. The subgraph in Figure 6a is the locally enlarged image of experimental ${\beta }_2$ by narrowing the scale range of vertical direction. We can see that ${\beta }_0$ is strong and linearly proportional to the propagation distance. It grows without the limitation of phase-velocity matching. On the other hand, ${\beta }_2$ is weak and has a maximum cumulative propagation distance ($L_n$) due to the $\pi$ effect caused by the slight deterioration of phase-velocity matching condition ($L_n=$ $2\pi /(k_n^{*}-2k)$ for $k_n^{*}\ne 2k$ [17,28]. And as shown in Figure 1d, the condition of phase-velocity matching is not strictly satisfied, when considering S0 (400 kHz) as the second harmonic of the primary mode. This phase-velocity mismatching ($k_n^{*}\ne 2k$) may cause a periodic fluctuation of the second harmonic. The similar experimental results of the second harmonic have been well documented in the literature [17]. ${\beta }_0$ of the zero-frequency mode is much stronger than ${\beta }_2$ of the second harmonic at all signal sampling points. To investigate the robustness of ${\beta }_0$, different frequencies of experimental excitation signals, for example, 250 and 300 kHz are considered and the similar results are obtained shown in Figure 7. $\beta$ parameters are also presented for other frequencies in literature [39] and the numerical results of literature [39] are consistent with our theoretical predictions. It is confirmed that the zero-frequency mode can be utilized for evaluating the early-stage material nonlinearity.

To verify that ${\beta }_0$ can be used to monitor local damage, a simulation model as shown in Figure 8 is built. The material properties of Parts I and III are shown in Table 1 and the material properties of Part II is shown in

Table 2. Material properties of Al in Part II.

ρ (kg/m3)	λ (GPa)	μ (GPa)	$\mathcal{A}\left(\mathbf{GPa}\right)$	$\mathcal{B}\left(\mathbf{GPa}\right)$	$\mathcal{C}\left(\mathbf{GPa}\right)$
2704	70.3	26.96	−41600	−13100	−15050

. It can be noted that the three-order elastic coefficients in Table 2 is much larger than the ones in Table 1. Thus, it can be considered that the model contains the local fatigue damage in Part II. The width of Part II is 200 mm. We consider two locations of Part II: 700 and 800 mm away from the left end of the model, to investigate the change of ${\beta }_0$. The same nineteen signal detection points with the ones in Section 4 are considered. Figure 9 shows the simulation results of ${\beta }_0$ and ${\beta }_2$ for comparison. It can be seen that both ${\beta }_0$ and ${\beta }_2$ have a jump when waves propagate through the region of Part II and ${\beta }_0$ is much more sensitive to local change of the material nonlinearity. Therefore, the zero-frequency mode can be used as an effective efficient evaluation index for the early-stage material nonlinearity.

6. Conclusions

In conclusion, we present theoretical analysis and experimental and simulation results to demonstrate that the symmetric zero-frequency mode is effective for evaluating the early-stage material nonlinearity. Compared with the conventional techniques based on the second or higher-order harmonics of Lamb waves, the zero-frequency mode has many advantages, such as strong signals, easy accumulation, long-range inspection, arbitrary excitation frequency and no requirement of phase-velocity matching.

Theoretically, besides S0, other wave modes, for example, S1, A0, A1, etc. (Figure 1d), can also be used as the primary wave and yet only the symmetric zero-frequency mode can be generated due to the requirement of the non-zero power flux. The above method of using the symmetric zero-frequency mode is also valid for the case of using broad-band signals as the primary waves, which can be conveniently and inexpensively excited. Besides the weak material nonlinearity, we have also identified that the zero-frequency mode can be induced by a micro crack in an elastic plate, leading to the possibility of using the zero-frequency mode to detect and evaluate the micro crack or numerous distributed micro cracks. The above issues need to be confirmed with more experimental evidences in future. Nevertheless, the results presented in this work support that the zero-frequency mode can be used as an effective efficient evaluation index for the early-stage material nonlinearity.