Second-Harmonic Generation in Membrane-Type Nonlinear Acoustic Metamaterials

Abstract

We study analytically and numerically the second-harmonic generation in a one-dimensional nonlinear acoustic metamaterial, composed of an air-filled waveguide periodically loaded by clamped elastic plates. Based on the transmission line approach, we derive a nonlinear dynamical lattice model which, in the continuum approximation, leads to a nonlinear dispersive wave equation. By applying the perturbation method to the latter, we derive the analytical expressions for the first- and second-harmonics, which are in excellent agreement with the numerical simulations of the nonlinear dynamical lattice model. Apart from the case of dispersionless nonlinear propagation and the Fubini solution, special attention is payed to the role of dispersion. In that regard, it is found that, once dispersion comes into play, second-harmonic beatings in space due to phase-mismatch can be identified. Our results provide many opportunities for the development of new periodic acoustic structures featuring both nonlinearity and dispersion.

1. Introduction

Phononic crystals and acoustic metamaterials have been widely popularized since the beginning of the 21st century because of the exotic properties they offer, including subwavelength focusing, cloaking, and extraordinary transmission among others [1]. These materials are composite structures designed to tailor acoustic wave dispersion through Bragg’s scattering and local resonances. One of the main properties of the phononic crystals and acoustic metamaterials is that—thanks to their periodic structure—they exhibit phononic band gaps, namely ranges of frequencies where no propagation occurs (i.e., linear acoustic waves are evanescent). The existence of band gaps has been studied in theory and observed experimentally for the first time in periodic acoustic waveguides by Sugimoto [2] and Bradley [3]. In 2000, Liu et al. [4] paved the way to acoustic metamaterials through phononic crystals that exhibited spectral gaps with lattice constants two orders of magnitude smaller than the relevant acoustic wavelength. The formation of band gaps in these acoustic metamaterials is based on the idea of the inclusion of locally resonant structures. These, in turn, usually determine the properties of acoustic metamaterials, rather than their composition.

Most of the works in the field of phononic crystals and acoustic metamaterials are restricted in the linear regime, and they do not consider the nonlinearity of the medium. Nevertheless, as the amplitude of the wave excitation is increased, the response of the metamaterial becomes nonlinear. This may give rise to different phenomena including, for instance, harmonic generation [5,6] and emergence of solitons [7], namely robust localized waves propagating undistortedly due to a balance between dispersion and nonlinearity. Nonlinear acoustic metamaterials are very good candidates to analyze the combined effects of nonlinearity and dispersion, occurring, e.g., in the beating of higher generated harmonics, because of mismatched phases or the existence of solitons [8,9].

In the context of electromagnetic (EM) metamaterials, there exist many works devoted to the nonlinear behavior [10,11,12,13,14]. Typically, metamaterials can be realized or modeled by a quasi-lumped transmission line (TL), with elementary cells consisting of a series inductor and a shunt capacitor, the dimensions of which are much less than the wavelength of the operating frequency. The TL approach is a powerful tool for studying nonlinear phenomena in EM metamaterials, such as soliton formation and nonlinear propagation [11,12,13,14]. In the context of acoustic metamaterials, one may similarly employ an acoustic circuit modeling, in which the voltage corresponds to the acoustic pressure and the current to the volume velocity flowing through the waveguide, in order to characterize the nonlinear propagation. On the other hand, the linear TL description of acoustic metamaterials has gained considerable attention the last few years [15,16,17,18,19]. However, studies on nonlinear phenomena in such settings are rather limited [9].

In this work, we consider an acoustic metamaterial composed of an air-filled waveguide, periodically side-loaded by clamped elastic plates. For such a system, it is well known that the elastic plates are incorporated as resonant elements and are considered to be in series in the electro-acoustic analogy [15,16,17,18,19]. Nonlinear wave propagation in this dispersive structure is theoretically and numerically studied. The paper is structured as follows. Based on the nonlinear transmission line theory, in Section 2, we introduce the 1D nonlinear lattice model, and employing the continuum approximation, we derive the nonlinear dispersive wave equation. Linear and nonlinear properties of the model, namely the dispersion relation and the case of dispersionless nonlinear propagation, are respectively studied in Section 3 and Section 4. In Section 5, by applying a perturbation method, we study the second-harmonic generation, while the effect of dispersion is studied in detail in Section 6. Finally, in Section 7, we present our conclusions and discuss future research directions.

2. Electro-Acoustic Analogue Modeling

We consider low-frequency wave propagation in an acoustic waveguide periodically loaded with elastic plates. The frequency range considered is well below the first cut-off frequency of the higher propagative modes in the waveguide, therefore the problem is considered as one-dimensional. The distance between the plates is d and the plates have a thickness h and radius r, as shown in Figure 1a.

In order to theoretically analyze this system, in this work, we adopt the electro-acoustic analogy. Using physical acoustics, one has to treat our setting by solving two nonlinear partial differential equations for the pressure and velocity field coupled at specific points (where the resonators are located) with a number of ordinary differential equations that describe the dynamics of the resonators (in our case the clamped elastic plates). This kind of modeling is very hard to treat analytically and one has only to rely on numerical simulations. On the contrary, one could use the electro-acoustical analogy to derive a nonlinear discrete wave equation, describing wave propagation in an equivalent electrical transmission line, which can be solved perturbatively in the continuum limit. Such an approach provides an efficient way to treat the nonlinearity and greatly simplifies the problem, allowing for straightforward analytical treatment by means of standard techniques that are used in other physical systems. Then, in this work, the voltage v and the current i of the equivalent electrical transmission line corresponds to the acoustic pressure p and to the volume velocity u flowing through the waveguide cross-section, respectively.

Following the TL approach, we start our consideration with the unit-cell circuit of the equivalent TL model of this setting. There are two different forms for the unit-cell circuit, as shown in Figure 1b and Figure A1 (see Appendix). Here, we only introduce the first one. The other circuit is introduced in the Appendix [18,19].

The unit-cell is composed of two parts, one corresponding to the propagation in the tube and the other one to the elastic plate. The resonant elastic plate can be modeled by an $LC$ circuit, namely the series combination of an inductance $L_m$ and a capacitance $C_m$, given by

$$L_m=\frac{{\rho }_{m}h}{S},C_m=\frac{1}{{\omega }_m^2{L}_m},$$

(1)

where ${\rho }_m$ is the plate density, S represents the cross-section area of the plate, while ${\omega }_m=2\pi f_m$ is the resonance frequency of the plate, with

$$f_m=0.4694\frac{h}{r^2}\sqrt{\frac{E}{{\rho }_m(1-{\nu }^2)}},$$

(2)

where r is the radius of the elastic plate, E is the Young’s modulus and ν is the Poisson ratio [18,19]. We consider elastic plates made of rubber, with ${\rho }_m=1420$ kg/m${}^3$, $E=2.758$ GPa and $\nu =0.34$. The part of the unit-cell circuit that corresponds to the waveguide is modeled by the inductance $L_{\omega }$ and shunt capacitance $C_{\omega }$; the linear parts of these elements are given by

$$L_{\omega 0}=\frac{{\rho }_{0}d}{S},C_{\omega 0}=\frac{Sd}{{\rho }_0{c}_0^2},$$

(3)

where ${\rho }_0$ and $c_0$ are, respectively, the density and the sound velocity of the fluid in the waveguide that has a cross section $S=\pi r^2$.

In this work, we consider the response of the elastic plate to be linear and the propagation in the waveguide weakly nonlinear. This is a reasonable approximation, since the amplitudes used in this work are not sufficient to excite nonlinear vibrations of the elastic plate [9,20]. Due to the compressibility of the air, the wave celerity is a nonlinear term, $c_{NL}$. Thus, we consider that the capacitance $C_{\omega }$ is nonlinear and depends on the pressure p. Approximating the celerity as $c_{NL}\approx c_0\left(1+{\beta }_{0}p/{\rho }_0{c}_0^2\right)$, where ${\beta }_0$ is the nonlinear parameter for the case of air (that is ${\beta }_0=1.2$), the pressure-dependent capacitance $C_{\omega }$ can be expressed as

$$C_{\omega }=C_{\omega 0}-C_{\omega }^{{}^{\prime }}{p}_n,$$

(4)

$$C_{\omega }^{{}^{\prime }}=\frac{2{\beta }_0}{{\rho }_0{c}_0^2}{C}_{\omega 0}.$$

(5)

The inductance remains in the same form as in the linear part, $L_{\omega 0}=L_{\omega }$. Next, we will derive an evolution equation for the pressure in the n-th cell of the lattice as follows.

First, we note that the advantage of the considered unit-cell circuit is that the inductances $L_{\omega }$ and $L_m$ are in a series connection and, thus, can be substituted by the inductance $L=L_{\omega }+L_m$ (see Figure 1b). Applying Kirchoff’s voltage law for two successive cells yields

$$p_{n-1}-p_n=L\frac{d}{dt}{u}_n+V_n,$$

(6)

$$p_n-p_{n+1}=L\frac{d}{dt}{u}_{n+1}+V_{n+1},$$

(7)

where $V_n$ is the voltage produced by the capacitance of the elastic plates $C_m$. Subtracting the two equations above, we obtain the differential-difference equation (DDE)

$${\widehat{\delta }}^2{p}_n=L\frac{d}{dt}\left(u_n-u_{n+1}\right)+\left(V_n-V_{n+1}\right),$$

(8)

where ${\widehat{\delta }}^2{p}_n\equiv p_{n+1}-2p_n+p_{n-1}$. Then, Kirchhoff’s current law yields

$$u_n-u_{n+1}=C_{\omega }\frac{d}{dt}\left(p_n\right),$$

(9)

with

$$u_n=C_m\frac{d}{dt}\left(V_n\right)\text{and}{u}_{n+1}=C_m\frac{d}{dt}\left(V_{n+1}\right).$$

(10)

Subtracting Equation (10) and employing Equation (9), we obtain

$$u_n-u_{n+1}=C_m\frac{d}{dt}\left(V_n-V_{n+1}\right)=C_{\omega }\frac{d}{dt}\left(p_n\right).$$

(11)

Then, recalling that the capacitance $C_{\omega }$ depends on the pressure (cf. Equation (4)), we express $V_n-V_{n+1}$ as

$$V_n-V_{n+1}=\frac{C_{\omega }}{C_m}{p}_n=\frac{C_{\omega 0}}{C_m}{p}_n-\frac{C_{\omega }^{{}^{\prime }}}{C_m}{p}_n^2.$$

(12)

Next, substituting Equations (9) and (12) into Equation (8), we obtain the following evolution equation for the pressure

$${\widehat{\delta }}^2{p}_n=L\frac{d}{dt}\left(C_{\omega }\frac{d}{dt}\left(p_n\right)\right)+\frac{C_{\omega 0}}{C_m}{p}_n-\frac{C_{\omega }^{{}^{\prime }}}{C_m}{p}_n^2.$$

(13)

To this end, employing Equation (4), we can rewrite the above equation as follows:

$${\widehat{\delta }}^2{p}_n-L{C}_{\omega 0}\frac{d^2{p}_n}{d{t}^2}-\frac{C_{\omega 0}}{C_m}{p}_n+\frac{L{C}_{\omega }^{{}^{\prime }}}{2}\frac{d^2{p}_n^2}{d{t}^2}+\frac{C_{\omega }^{{}^{\prime }}}{C_m}{p}_n^2=0.$$

(14)

In this article, we have numerically integrated the nonlinear lattice model, Equation (14), by using the function ode45 of Matlab which is based on the Runge Kutta method, with an initial condition, $p_1=sin\left(\omega t\right)$ at $x=0$, where $\omega =2\pi f$ is the angular frequency of the driver. For each simulation, we ensure the validity of the Courant–Friedrichs–Lewy (CFL) condition, $c\frac{dt}{dx}\le 1$, where c is the phase velocity, $dt$ and $dx$ are the time step and length interval, respectively. We also pay attention to the length of the system, which should be long enough to avoid reflections in the analyzed signal.

2.1. The Continuum Approximation

For our analytical considerations, we will focus on the continuum limit of Equation (14), corresponding to $n\to \infty$ and $d\to 0$ (but with $nd$ being finite); in such a case, the pressure becomes $p_n\left(t\right)\to p(x,t)$, where $x=nd$ is a continuous variable, and

$$p_{n\pm 1}=p\pm d\frac{\partial p}{\partial x}+\frac{d^2}{2}\frac{{\partial }^{2}p}{\partial x^2}\pm \frac{d^3}{3!}\frac{{\partial }^{3}p}{\partial x^3}+\frac{d^4}{4!}\frac{{\partial }^{4}p}{\partial x^4}+O\left(d^5\right),$$

(15)

i.e., the difference operator ${\widehat{\delta }}^2$ is approximated by ${\widehat{\delta }}^2{p}_n\approx d^2{p}_{xx}+\frac{d^4}{12}{p}_{xxxx}$. Keeping the $O\left(d^4\right)$ derivative term, the PDE contains also the weak dispersion that originates from the periodicity of the elastic plates array as we will see later on. Terms of the order $O\left(d^5\right)$ and higher are neglected, and subscripts denote partial derivatives. This way, Equation (14), becomes the following PDE:

$$d^2{p}_{xx}+\frac{d^4}{12}{p}_{xxxx}-L{C}_{\omega 0}{p}_{tt}-\frac{C_{\omega 0}}{C_m}p+\frac{1}{2}L{C}_{\omega }^{{}^{\prime }}{\left(p^2\right)}_{tt}+\frac{C_{\omega }^{{}^{\prime }}}{C_m}{p}^2=0.$$

(16)

It is also convenient to express our model in dimensionless form; this can be done upon introducing the normalized variables τ and χ and normalized pressure p, which are defined as follows: τ is time in units of ${\omega }_B^{-1}$, where ${\omega }_B=\pi c_0/d$ is the Bragg frequency; χ is space in units of $c/{\omega }_B$, where the velocity is given by

$$c=\frac{c_0}{\sqrt{1+\alpha }},\alpha =\frac{h{\rho }_m}{d{\rho }_0},$$

(17)

and p is pressure in units of $p_0={\rho }_0{c}_0^2$. Then, Equation (16) is reduced to the following dimensionless form:

$$p_{\tau \tau }-p_{\chi \chi }-\gamma p_{\chi \chi \chi \chi }+m^{2}p=2{\beta }_0{m}^2{p}^2+{\beta }_0{\left(p^2\right)}_{\tau \tau },$$

(18)

where parameters $m^2$ and γ are given by

$$m^2=\frac{\alpha }{1+\alpha }{\left(\frac{{\omega }_m}{{\omega }_B}\right)}^2,\gamma =\frac{1}{12}{\pi }^2(1+\alpha ).$$

(19)

It is interesting to identify various limiting cases of Equation (18). First, in the linear limit $({\beta }_0=0$, or $p^2\ll 1)$, and in the absence of plates ($m^2\to 0$, and without considering higher order spatial derivatives), Equation (18) is reduced to the linear wave equation, $p_{\tau \tau }-p_{\chi \chi }=0$. In the linear limit, in the presence of plates and in the long wavelength approximation ($k\to 0$, and without considering higher order spatial derivatives), Equation (18) takes the form of the linear Klein–Gordon equation [7], $p_{\tau \tau }-p_{\chi \chi }+m^{2}p=0$, with the parameter m playing the role of mass. Finally, in the nonlinear regime, but when plates are absent, Equation (18) is reduced to the well-known Westervelt equation, $p_{\tau \tau }-p_{\chi \chi }-{\beta }_0{\left(p^2\right)}_{\tau \tau }=0$, which is a common nonlinear model describing 1D acoustic wave propagation [21].

3. The Linear Dispersion Relation

We now consider the linear limit of Equation (18) and the respective dispersion relation. Assuming propagation of plane waves, of the form $p\propto exp\left[i\right(k\chi -\omega \tau \left)\right]$, we obtain the following dispersion relation connecting the wavenumber k and frequency ω

$$D(\omega ,k)=-{\omega }^2+k^2-\gamma k^4+m^2=0.$$

(20)

For $D(\omega ,k)=0$, this is the familiar dispersion relation of the linear Klein–Gordon model. It is clear that Equation (20) suggests the existence of a gap at low frequencies, i.e., for $0\le \omega <m$, with the cut-off frequency defined by the parameter m. For $m<\omega <{\omega }_B$, there exists a band, with the dispersion curve $\omega \left(k\right)$ having the form of hyperbola, which asymptotes (according to Equation (20)) to unity, which is the normalized velocity associated with the wave equation $p_{\tau \tau }-p_{\chi \chi }=0$. The term $\gamma k^4$ accounts for the influence of the periodicity of the lattice (originating from the term ${\widehat{\delta }}^2{p}_n$) to the dispersion relation. Although this term appears to lead to instabilities for large values of k, both Equations (18) and (20) are used in our analysis only in the long wavelength limit where k is sufficiently small.

Since all quantities in the above dispersion relation are dimensionless, it is also relevant to express it in physical units. In particular, taking into account that the frequency ${\omega }_{ph}$ and wavenumber $k_{ph}$ in physical units are connected with their dimensionless counterparts through $\omega ={\omega }_{ph}/{\omega }_B$ and $k=\frac{k_{ph}c}{{\omega }_B}$, we can express Equation (20) in the following form:

$$-{\omega }_{ph}^2+k_{ph}^2{c}^2-\gamma \frac{k_{ph}^4{c}^4}{{\omega }_B^2}+m^2{\omega }_B^2=0.$$

(21)

Solving Equation (21) analytically with respect to $k_{ph}$, we can then determine the frequency $f={\omega }_{ph}/2\pi$ as a function of the wavenumber $k_{ph}$, and plot the resulting dispersion relation. The real and imaginary parts of the dispersion relation are respectively plotted in Figure 2a,b for a metamaterial composed of elastic plates made of rubber with thickness $h=2.78\times {10}^{-4}$ m and with a periodicity $d=0.01$ m. The dispersion relation features the band gap from 0 Hz to $\left(m\frac{{\omega }_B}{2\pi }\right)$ Hz due to the combined effect of the resonance of the plate and of the geometry of the system. The upper limit of this band gap is found to be sufficiently smaller than the Bragg band frequency $f_B=c_0/2d=17163$ Hz, with $c_0=343.26$ m/s. The propagating band has two parts: a strongly dispersive and a weakly dispersive one. In the lower weakly dispersive region, there is a “quasi-linear” dispersion with the slope $\mathrm{a}=\frac{c_0}{\sqrt{1+\alpha }}$ (which is identical to the velocity c in Equation (17)), and the upper weakly dispersive region is due to the periodicity of lattice. Both the periodicity of the system d and the thickness of the elastic plates could influence the first cut-off frequency m and the slope $\mathrm{a}$ of the “quasi-linear” dispersion, as shown in the Figure 2c,d. The first cut-off frequency m is inversely proportional to the periodicity of the lattice d and proportional to the thickness of the elastic plates h, while the slope $\mathrm{a}$ of the “quasi-linear” dispersion increases with the increase of d and the decrease of h. Due to periodicity, the band structure of our system exhibits a Bragg band gap with an upper edge $kd=\pi$ located at $17.163$ kHz. The lower edge of the gap, however, also depends on α (describing the impedance mismatch and the filling fraction) and is located much lower at $1.988$ kHz. Due to the dispersion around this lower band gap edge, the 4th order spatial derivative term is needed to describe the system in a good accuracy. To further illustrate the importance of the higher order dispersive term, in Figure 2a we additionally show a curve corresponding to the case without it ($\gamma =0$).

On the other hand, the red lines in the Figure 2a,b show the respective results (for the lossless case under consideration) for the dispersion relation, as obtained using the transfer matrix method (TMM) [3]

$$cos\left(k_{ph}d\right)=cos(\frac{{\omega }_{ph}}{c_0}d)+i\frac{Z_m}{2Z_0}sin(\frac{{\omega }_{ph}}{c_0}d),$$

(22)

where $Z_m=i\left({\omega }_{ph}{L}_m-1/{\omega }_{ph}{C}_m\right)$ is the impedance of the plate for the lossless case in the long wavelength approximation, and $Z_0={\rho }_0{c}_0/S$ the acoustic characteristic impedance of the waveguide. Comparing the dispersion relation obtained by using TMM, with the one resulting from the continuum approximation, Figure 2a,b, we find an excellent agreement between these two in the regime of low frequencies, namely sufficiently lower to the Bragg frequency.

4. Dispersionless Nonlinear Propagation

We now consider the nonlinear regime, without dispersion due to the periodicity of the lattice and the resonance of the plates, i.e., the well-known Westervelt equation. As the large-amplitude wave propagates, the amplitude of the fundamental component $|p_1|$ will decrease continuously as the energy is transferred to the nonlinearly generated higher-harmonic components ($|p_2|$, etc. ). The growth of the higher harmonics is displayed in the pre-shock region which is defined by $\sigma \le 1$, where $\sigma =x/x_{sh}$ is a dimensionless shock formation distance. The shock distance,

$$x_{sh}=\frac{\rho c^3}{2\pi f\beta p_0},$$

(23)

is proportional to the velocity and inversely proportional to the pressure amplitude and source frequency for a fixed medium. Here, the source condition is $p(0,t)=p_{0}sin\omega t$, where $p_0=0.04P_0$ (with $P_0$ being the atmospheric pressure), and $sin\omega t$ is periodic in time with fundamental frequency $f=400$ Hz. If this initial state propagates in a dispersionless waveguide—cf. Figure 3—its shock distance will be around 4 m. In the near source region, the pertinent Fubini solution has certain asymptotic properties [21]. The amplitude of the second harmonic component increases linearly to the propagation distance:

$$|p_1|/|p_0|=1-\frac{1}{8}{\sigma }^2+O\left({\sigma }^4\right),$$

(24)

$$|p_2|/|p_0|=\frac{1}{2}\sigma +O\left({\sigma }^3\right).$$

(25)

As shown in Figure 3, the numerical results (circles) are in a good agreement with the Fubini solution (solid lines).

5. Combining Dispersion and Nonlinearity: Perturbation Method

Now, we study the second-harmonic generation in the presence of the periodic array of the elastic plates, namely in the presence of dispersion. Our analysis relies on the determination of approximate solutions of Equation (18) by using a perturbative approach.

We now express p as an asymptotic series in ϵ

$$p=ϵp_1+{ϵ}^2{p}_2+{ϵ}^3{p}_3+\dots ,$$

(26)

where $0<ϵ\ll 1$ is a formal small parameter. Here, the introduced ϵ is the acoustic Mach number, defined as $ϵ=\frac{p_0}{\rho c^2}$, where $p_0$ is the amplitude of the incident wave. Then, substituting Equation (26) into Equation (18), we obtain a hierarchy of equations at various orders in ϵ. Of particular importance in our analysis are the equations at the first two orders, which are as follows.

At the leading order, $O\left({ϵ}^1\right)$, the resulting equation is

$$\frac{{\partial }^2{p}_1}{\partial {\tau }^2}-\frac{{\partial }^2{p}_1}{\partial {\chi }^2}-\gamma \frac{{\partial }^4{p}_1}{\partial {\chi }^4}+m^2{p}_1=0,$$

(27)

which possesses plane wave solution of the form

$$p_1=\frac{1}{2}exp\left(i\theta \right)+\mathrm{c}.\mathrm{c}.\equiv cos\left(\theta \right),$$

(28)

where c.c. denotes complex conjugate, $\theta =\omega \tau -k\left(\omega \right)\chi$ is the phase, while parameters k and ω satisfy the dispersion relation Equation (20). Next, we consider the equation at the order $O\left({ϵ}^2\right)$,

$$\frac{{\partial }^2{p}_2}{\partial {\tau }^2}-\frac{{\partial }^2{p}_2}{\partial {\chi }^2}-\gamma \frac{{\partial }^4{p}_2}{\partial {\chi }^4}+m^2{p}_2=2{\beta }_0{m}^2{p}_1^2+{\beta }_0\frac{{\partial }^2{p}_1^2}{\partial {\tau }^2}.$$

(29)

Substituting Equation (28) into Equation (29), and using the identity ${cos}^2\left(\theta \right)=\left(cos\left(2\theta \right)+1\right)/2$, we rewrite Equation (29) as follows:

$$\frac{{\partial }^2{p}_2}{\partial {\tau }^2}-\frac{{\partial }^2{p}_2}{\partial {\chi }^2}-\gamma \frac{{\partial }^4{p}_2}{\partial {\chi }^4}+m^2{p}_2=\left(m^2-2{\omega }^2\right){\beta }_{0}cos\left(2\theta \right).$$

(30)

The solution of this equation is the sum of the general solution $p_2^h$ of the homogeneous equation and the particular solution $p_2^p$ of the inhomogeneous equation, namely $p_2=p_2^h+p_2^p$, where the corresponding waves for these two solutions are the free and forced waves respectively; these solutions read

$$\begin{array}{c}\hfill p_2^h=p_2^h(x=0)cos\left(2\varphi \right)=\frac{1}{2}{p}_2^h(x=0)exp\left(i2\varphi \right)+\mathrm{c}.\mathrm{c}.\end{array},$$

(31)

$$\begin{array}{c}\hfill p_2^p=\frac{m^2-2{\omega }^2}{D(2\omega ,2k)}{\beta }_{0}cos\left(2\theta \right)=\frac{1}{2}\frac{m^2-2{\omega }^2}{D(2\omega ,2k)}{\beta }_{0}exp\left(i2\theta \right)+\mathrm{c}.\mathrm{c}.\end{array},$$

(32)

where $2\varphi =2\omega \tau -k\left(2\omega \right)\chi$, and $k\left(2\omega \right)$ is the wavenumber of the free wave at second harmonic frequency. As long as $2k\left(\omega \right)\ne k\left(2\omega \right)$, which is the case in dispersive media, the forced and free waves have different phase speeds, i.e., they are phase-mismatched. Since there is no second harmonic at $x=0$, we can set

$$p_2^h(x=0)=-\frac{m^2-2{\omega }^2}{D(2\omega ,2k)}{\beta }_0.$$

(33)

Thus, the evolution of the second harmonic field $p_2$ can directly be found as a combination of Equations (32) and (31):

$$p_2=-2\frac{m^2-2{\omega }^2}{D(2\omega ,2k)}{\beta }_{0}sin\left(\frac{\Delta k}{2}\chi \right)sin(2\omega \tau -k_{\mathrm{eff}}\chi ),$$

(34)

where $k_{\mathrm{eff}}$ is the effective wave number,

$$k_{\mathrm{eff}}=\left(2k\left(\omega \right)+k\left(2\omega \right)\right)/2,$$

(35)

while $\Delta k$ is the detuning parameter that describes the asynchronous second harmonic generation,

$$\Delta k=k\left(2\omega \right)-2k\left(\omega \right)=k_2-2k_1.$$

(36)

Obviously, in the linear limit (${\beta }_0=0$), $p_2$ turns to 0, i.e., the generated second harmonic is due to the nonlinear effect. We can also find the second harmonic beatings in space, $sin\left(\frac{\Delta k}{2}\chi \right)$ in Equation (34). The position of the maximum of the beating can be related to the second-harmonic phase-mismatching frequency as

$$x_c\left(n\right)=\frac{\pi }{\Delta k_n}=\frac{\pi }{\left|k\left(n\omega \right)-nk\left(\omega \right)\right|}.$$

(37)

Therefore, as $\Delta k$ increases, the second harmonic beating spatial period, and also its maximum amplitude, decreases.

6. Results

Now, we study numerically the role of dispersion on harmonic generation, in the 1D acoustic metamaterial composed of elastic plates, and compare the numerical results to the analytical findings of Section 5. There are two cases, corresponding to the propagating driver and evanescent driver, which will be investigated separately below.

6.1. Driving Frequency in the Pass Band

We start by studying the case where both the fundamental component and the second harmonic component are in the pass band. We numerically integrate the weakly nonlinear lattice model, Equation (14), using an initial condition $p_1(x=0)$ of a cosinusoidal form, and also determine the spectrum of the solution by using the Fast Fourier Transform (FFT). Blue circles and red squares in Figure 4 show the evolution of the amplitude of the fundamental and second-harmonic components, respectively, as the wave propagates in the dispersive structure, as obtained numerically. As shown in the last section, if the distance between two cells is very small (here, we choose $d=0.01$ m), the asymptote, $\mathrm{a}$, in such a dispersive system will be very small—cf. green dotted line in Figure 2c—thus, the shock distance has a very small value. This means that the higher-harmonic generation process is achieved much closer to the source, compared to the Fubini solution (see Figure 3). The second-harmonic component $\left|p_2\right|/\left|p_0\right|$ no longer increases linearly because of the presence of dispersion induced by the elastic plates. Notice that $\left|p_2\right|/\left|p_0\right|$ develops beatings in space due to the phase mismatch. The fundamental wave vector of our source $k_1$ (frequency ω) generates a forced wave $2k_1$ (frequency $2\omega$), while the free wave that our system allows to propagate is $k_2$. The difference between $2k_1$ and $k_2$ (because of the dispersion relation) introduces a phase mismatch, which means that the second harmonic generation is asynchronous. We can, therefore, clearly observe the beatings in Figure 4, where the analytical solutions have an excellent agreement with the numerical ones. This perfect agreement breaks down if we consider a driving frequency with a second harmonic higher than 1 kHz, since the continuous dispersion relation, see Equation (20), deviates from the discrete one. Additionally, since we have considered a weakly nonlinear regime, the agreement between numerics and our approximation is found to break down for amplitudes larger than 10 kPa. When we increase the frequency from 350 Hz to 400 Hz, $\Delta k$ decreases, $x_c$ increases, the second-harmonic beatings spatial period—and also its maximum amplitude—increase, see Equation (34). During the nonlinear propagation in the dispersive system, cumulative nonlinear effects generate harmonics of the fundamental frequency, and we can control this process by tuning the dispersion relation with either the properties of the array or that of the plates.

6.2. Driving Frequency in the Bandgap

When the driving frequency is in the band gap, its second harmonic may be located either in the gap band (evanescent) or in the pass band (propagating). In the former case, the second harmonic is damped, and its decay rate is given by the imaginary part of the dispersion relation. In the latter case, the second harmonic is propagating through the structure. Both cases are analyzed below.

We start by studying the case where both the fundamental component and its second harmonic are in the band gap. In this case, the dispersion relation does not support real solutions, so the corresponding wave number is imaginary, namely $k\left(\omega \right)=i{k}_1^{\prime \prime }$, with $k_1^{\prime \prime }$ the imaginary part of the wavenumber, given by the dispersion relation. When $k_1^{\prime \prime }$ is very large, for example at $f=90$ Hz, cf. Figure 5a, the fundamental component $\left|p_1\right|$ decreases exponentially very quickly, namely $\left|p_1\right|\propto exp[-\mathrm{Im}\left(k\right)x]$. The second harmonic component is generated at the beginning of the structure but is very small. In this case, $k\left(2\omega \right)=i{k}_2^{\prime \prime }$, with $k_2^{\prime \prime }$ the imaginary part of the wavenumber given by the dispersion relation for $2\omega$. The generated frequency $2f=180$ Hz is still in the band gap, the corresponding $\mathrm{Im}\left(k\right)$ is smaller than that of the fundamental frequency, but still very large. Therefore, we can find a very small $\left|p_2\right|/\left|p_0\right|$ at the beginning of the structure, which eventually decreases to zero, see Figure 5b.

If the driving frequency is close to the first cut-off frequency, the frequency of the generated second-harmonic will be in the pass band. In this case, $\left|p_1\right|/\left|p_0\right|$ decays slowly during the propagation because the imaginary part of k is very small, cf. Figure 5c, where $f=300$ Hz, while the generated second-harmonic component is propagating, cf. Figure 5d. We can also find the second-harmonic beatings at the beginning due to the phase mismatch (the difference between $2k_1$ and $k_2$). Comparing to the result shown in Figure 4, the beating spatial period is smaller because of a bigger $\Delta k$. The fundamental component is evanescent and, therefore, there will not be phase mismatch for $\left|p_2\right|$ (only $k_2$ propagates), i.e., we can hardly observe the beatings for $\left|p_2\right|/\left|p_0\right|$ after $0.4$ m, as shown in Figure 5d.

7. Conclusions

In conclusion, we have theoretically and numerically studied the nonlinear propagation and second-harmonic generation in 1D acoustic metamaterial composed of an air-filled tube with a periodic array of elastic plates. Based on the electro-acoustic analogy, we proposed the transmission line approach to derive a nonlinear lattice model, which was analyzed by both numerical and analytical techniques. Considering the continuum limit of the lattice model, we derived a nonlinear dispersive wave equation, in the form of a nonlinear Klein–Gordon model, which reduces—at certain limits—to other well-known acoustic wave models (such as the Westervelt equation). In the linear limit, we derived from this model the dispersion relation which, in the low frequency regime, was found to be in excellent agreement with the one obtained by the transfer matrix method. We have shown that, during the nonlinear propagation, cumulative nonlinear effects generate harmonics of the fundamental frequency. Dispersion introduces phase mismatch between at higher harmonics, which is the responsible of the beating effect. We used a perturbative approach to study analytically the effect of dispersion on the harmonic generation. Analytical and numerical results were found to be in excellent agreement.

There are many future research directions that may follow these preliminary results. First, it would be interesting to investigate if the combined effects of nonlinearity and dispersion may give rise to the emergence of solitons in the system (i.e., nonlinear waves that propagate undistorted when nonlinearity and dispersion are exactly counterbalanced). Second, taking into account the presence of inherent losses in the metamaterial structures under consideration, a study of the interplay between nonlinearity, dispersion and losses would be very relevant. It would also be interesting to extend our work to higher-dimensional settings, as well as to study nonlinear propagation in double negative metamaterials, waveguides periodically loaded with a combination of elastic plates and Helmholtz resonators.