Elastic wave propagation in nonlinear two-dimensional acoustic metamaterials

This study describes the wave propagation in a periodic lattice which is formed by a spring-mass two-dimensional structure with local Duffing nonlinear resonators. The wave propagation characteristics of the system are evaluated by using the perturbation method to determine the dispersion relationships and wave propagation characteristics in the nonlinear two-dimensional acoustic metamaterials. A quantitative study of wave amplitude is carried out to determine the maximum allowable wave amplitude for the whole structures under the assumption of small parameters. In particular, the harmonic balance method is introduced to investigate the frequency response and effective mass of the nonlinear systems. We find that the dispersion relations and group velocity of unit cell are related to wave amplitude. Furthermore, the dual-wave vector is observed in the nonlinear systems. Numerical simulations validate the dispersion analytical results. The results can be used to tune wave propagation in the nonlinear acoustic metamaterials and provide some ideas for the study of nonlinear metamaterials.


Introduction
Lightweight periodic lattices [1] have been developed and used in disparate engineering applications for many decades [2] owing to their distinguished properties such as high specific strength, high stiffness and multi-functional potential [3,4]. In addition to their superior static properties, lattice structures also provide excellent dynamic properties. The wave propagation behaviors in periodic media are interesting and important dynamic phenomena and properties. Owing to the periodic properties, periodic structures show extraordinary behaviors for prohibiting wave propagation, at locations called ''stop bands'' in a specific frequency range. After past decades of researches, two types of stop bands can be identified in periodic systems. One type is represented by Bragg band gaps, which are mainly related to the material properties or the geometry of the unit cell using properties like hierarchical structures [5], auxetics [6] and structures with curved beams [7]. Boundary conditions periodicity also contributes to Bragg band gaps [8]. The other types of stop bands are the socalled locally resonant (LR) band gaps, caused by a resonating system included in the periodic structure. The pioneering work of Liu et al. [9] has led to the generation of acoustic metamaterials (AM) which are typically artificial periodic media that has subwavelength characteristics of elastic waves based on mechanism of local resonance. However, band gap behavior is not limited to periodic structures, and it can also be observed in aperiodic structures [10]. At present, numerous studies focus on linear AM [11][12][13][14]. There has been an interest in designing AM to present negative mass effect [15]. Furthermore, LR band gaps can happen at wavelengths which are much larger than the lattice size, yielding low-frequency vibration/noise attenuation and wave filtering [16,17]. However, in such linear configurations, the resulting attenuation bandwidth and wave propagation in AM are limited by the parameters of material and structure, i.e., masses of local resonators [18], which is typically aimed to be minimized in most applications that require lightweight and the wave propagation of linear systems cannot be tuned in active or passive ways.
Compared with linear systems, nonlinearity provides abundant dynamic characteristics, such as subharmonics, solitary waves, wave coupling, phase matching, self-collimation, multi-stable, nonreciprocal diodes [19][20][21][22][23]. In addition, nonlinear dynamic phenomena such as nonlinear resonance, modal coupling, chaos [24][25][26] are expected to provide new wave control mechanisms for nonlinear acoustic/ elastic media. The introduction of nonlinear factors into AM is expected to break the technical bottleneck of low-frequency, broadband and high-efficiency vibration control of structures, which gradually become an emerging hot topic in the metamaterial research community.
At present, a few researchers [27][28][29] have reported the dispersion characteristics of nonlinear AM. In the field of energy harvesting, appropriately designed nonlinear oscillator can be used to enhance frequency bandwidth [30]. And a plethora of wideband dynamic behavior is discovered in bistable configuration that chaotic interwell vibrations depend on the input amplitudes [31,32]. As for the nonlinear vibration suppression, Higashiyama et al. [33] studied the bandedge mode of the diatomic chain model under evenorder nonlinearity. Cveticanin et al. [34] studied the equivalent negative mass characteristics of the Duffing nonlinear diatomic chain metamaterial model based on theory. Zhou et al. [35] established a diatomic chain metamaterial model with a linear oscillator attached to a nonlinear matrix structure. Ganesh [36] and Bukhari [37] worked on the spectrospatial analyses of nonlinear metamaterial to elucidate the relation between topological/physical (space-time domain) and spectral (wave number-frequency domain) features of wave motion. The space wave packet propagation characteristics in an infinite weakly nonlinear structure were studied, and it was found that the transient wave packet has a significant frequency shift when propagating in the band gap due to the nonlinearity. Gantzounis et al. [38] attached a local resonator to a sphere of a granular crystal to form a nonlinear metamaterial and studied its band gap and damping characteristics. Khajehtourian and Hussein [39] studied the approximate dispersion characteristics of a one-dimensional nonlinear locally resonant metamaterial rod based on the transfer matrix method. The nonlinearity of the system originates from the high-order deformation of the thin-walled rod. The position and width of the band gap vary with the wave amplitude. Bilal et al. [40] used the coupling of vibrator and magnetic force to induce nonlinearity and realized a transistor-like phonon switching element at a single frequency of one-dimensional metamaterial.
However, previous researches on nonlinear acoustic metamaterials mainly focused on one-dimensional system and rarely involved nonlinear two-dimensional acoustic metamaterials (2D AM). Fronk et al. [41] presented higher-order multiple scales analysis aimed at revealing angle-and amplitude-dependent invariant waveforms and plane-wave stability, in two-dimensional periodic media. Manktelow et al. [42] using a finite-element discretization of a single unit cell followed by a perturbation analysis studied the wave propagation in continuous, periodic structures subject to weak nonlinearities.
The nonlinear two-dimensional systems have numerous coupling degrees of freedom which makes it difficult to establish and solve the equation of motion, especially in a high-dimensional complex nonlinear system. This leads to a great gap in the study of nonlinear 2D AM. As well known, vibration suppression performance, static performance and wave propagation characteristics of AMs play an important role in applications. It is necessary to investigate these characteristics of nonlinear 2D AM. This paper focuses on the wave propagation, vibration and static characteristics of nonlinear 2D AM and organizes as follows. In Sect. 2, the Floquet-Bloch theorem and perturbation method are introduced to investigate the dispersion relationships for linear and nonlinear 2D AM. The peculiar phenomena accompanied with the amplitude-dependent iso-frequency contours and group velocity are also studied in Sect. 2. Section 3 mainly introduces the vibration response of nonlinear systems calculated by the harmonic balance method. The excitation amplitudedependent is investigated in the vibration response. The effective mass of the nonlinear 2D AM is studied in Sect. 4. Summarization and conclusion for nonlinear effects on the nonlinear 2D AM are given in Sect. 5.

Dispersion analysis of two-dimensional metamaterials
In this section, the calculation method of the dispersion relation of infinite two-dimensional metamaterials is mainly introduced. Firstly, the case of a linear spring-mass resonator is considered; then, further investigation is dedicated to the case of a nonlinear spring-mass resonator.

Linear spring-mass resonator
Let us consider an infinite 2D periodic AM with a lattice constant l consisting of AM unit cells and a series of springs as shown in Fig. 1a. The unit cells in our AM are arranged in a square configuration with each host medium connected to eight neighboring unit cells through linear stiffness springs, among which four are diagonal. As shown in Fig. 1b, the linear unit cell is composed of the resonator m 3 , the host medium m 1 and connected by a linear stiffness k 3 spring. As for the connection between the host mediums, it is assumed the stiffness of springs in the x-axis and yaxis directions as k 1 and that of the diagonal ones as k 2 .
The AM is infinite and the (n; p)th unit cell is defined to apply the Floquet-Bloch theorem. For the plane time harmonic propagation in the xaxis and y-axis directions with small displacements and the coupling of horizontal and vertical motion, the motion equation of the (n; p)th unit cell in the system is given by: m 1 o 2 u 1ðn;pÞ ot 2 ¼k 1 u 1ðnþ1;pÞ þ u 1ðnÀ1;pÞ À 2u 1ðn;pÞ À Á þ 1 2 k 2 u 1ðnþ1;pþ1Þ þ u 1ðnÀ1;pÀ1Þ À 2u 1ðn;pÞ þ v 1ðnþ1;pþ1Þ þ v 1ðnÀ1;pÀ1Þ À 2v 1ðn;pÞ À Á þ 1 2 k 2 u 1ðnÀ1;pþ1Þ þ u 1ðnþ1;pÀ1Þ À 2u 1ðn;pÞ À v 1ðnÀ1;pþ1Þ À v 1ðnþ1;pÀ1Þ þ 2v 1ðn;pÞ À Á þ k 3 u 3ðn;pÞ À u 1ðn;pÞ À Á ; ð1Þ m 3 o 2 u 3ðn;pÞ ot 2 ¼ k 3 u 1ðn;pÞ À u 3ðn;pÞ þ v 1ðnÀ1;pÀ1Þ À 2v 1ðn;pÞ þ u 1ðnþ1;pþ1Þ þ u 1ðnÀ1;pÀ1Þ À 2u 1ðn;pÞ À Á þ 1 2 k 2 v 1ðnÀ1;pþ1Þ þ v 1ðnþ1;pÀ1Þ À 2v 1ðn;pÞ À u 1ðnÀ1;pþ1Þ À u 1ðnþ1;pÀ1Þ þ 2u 1ðn;pÞ À Á þ k 3 v 3ðn;pÞ À v 1ðn;pÞ À Á ; ð3Þ m 3 o 2 v 3ðn;pÞ ot 2 ¼ k 3 v 1ðn;pÞ À v 3ðn;pÞ where n and p are arbitrary integer numbers, u 1 and v 1 are the displacements of the host medium in x-axis and y-axis directions, u 3 and v 3 are the displacements of the resonator in x-axis and y-axis directions, respectively. The wave solution is assumed as: u rðn;pÞ ¼ A r e iðnk x þpk y ÀxtÞ ðr ¼ 1; 3Þ; ð5Þ v rðn;pÞ ¼ B r e iðnk x þpk y ÀxtÞ ðr ¼ 1; 3Þ; ð6Þ where A r and B r are the wave amplitudes, i ¼ ffiffiffiffiffiffi ffi À1 p is the imaginary number, x is the wave frequency, and [k x ,k y ] is the two components of the wavevector k. Substituting the solution of Eqs. (5) and (6) into Eqs. (1)-(4), we get: Equations (7)-(10) form an eigenvalue problem governed by: where U ¼ ½A 1 ; A 3 ; B 1 ; B 3 T , Sðk x ; k y Þ is a 4 9 4 matrix and I is an identity matrix. By using the condition of nontrivial solutions for Eq. (11), the relation between the wave frequency x and wavevector components k x and k y can be obtained. Figure 2a shows the first Brillouin zone of 2D periodic structures, the wavevector should be taken in the whole Brillouin zone in usual, while the unit cell is square and symmetrical, and the analysis can be restricted to the irreducible Brillouin zones (shadow areas). Furthermore, it is only necessary to search the boundary of shadow areas, i.e., along the path M À C À X À M, since the extremums of the wave frequencies are always found on the zone boundary. The dispersion curve of linear AM is shown in Fig. 2b, with the specific parameters for We can find that a local resonance band gap exists between normalized frequency X ¼ 0:95 and X ¼ 1:4, including the natural frequency of the linear resonator (x 3 ¼ ffiffiffiffiffiffiffiffiffiffiffiffi k 3 =m 3 p ¼ 1) within its range.
As for the treatment of equations of the e 1 order, we need to find the relationship between u ð1Þ 3ðn;pÞ , v ð1Þ 3ðn;pÞ and u ð1Þ 1ðn;pÞ , v ð1Þ 1ðn;pÞ . In this way, we can obtain all the nonlinear terms in the same equation. For this purpose, the solution form of u ð1Þ 3ðn;pÞ and v ð1Þ 3ðn;pÞ can be expressed under the form: where D and E are the corresponding wave amplitudes, while they could have no effect on the following calculations. Substituting Eqs. (25) and (26) into Eq. (24), we can get: 1ðn;pÞ 1ðn;pÞ 1ðn;pÞ 3 3ðn;pÞ 1ðn;pÞ 3ðn;pÞ 1ðn;pÞ 1ðn;pÞ 1ðn;pÞ 1ðn;pÞ 3 3ðn;pÞ 1ðn;pÞ 3ðn;pÞ 1ðn;pÞ 1ðn;pÞ ¼ c 1 e iðnk x þpk y ÀsÞ þ c 2 e 3iðnk x þpk y ÀsÞ þ c:c: 1ðn;pÞ 1ðn;pÞ ¼ c 3 e iðnk x þpk y ÀsÞ þ c 4 e 3iðnk x þpk y ÀsÞ þ c:c: The terms c 2 andc 4 are associated with the e 3iðnk x þpk y ÀsÞ terms. The linear kernel of Eq. (27) is similar to the e 0 -equations in Eq. (23). This implies that terms in the right of Eq. (27), which present as the function of e iðnk x þpk y ÀsÞ , make to force the updated e 1equations at resonance, thus resulting in unbounded solutions. These terms are defined as secular terms [44], which should be eliminated to obtain a bounded solution at e 1 order. Therefore, the values of c 1 and c 3 must be equal to zero and we can get an expression for the frequency correction x 1 by assuming the wave amplitude A 1 ¼B 1 ,: The amplitude-dependent correction of the wave frequency has been obtained by many authors in different method, including the homotopy method [45], the multiple scale method [46,47] and the L-P method [48]. According to Eqs. (11), (22) and (29), each frequency X j is given as a summation of the linear zero order and perturbed system solution, as follows: Equation (30) gives the dispersion relationship of 2D AM with nonlinear resonator. It is indicated that the frequency is mainly affected by nonlinear stiffness k n , wave amplitude A 1 . For the following analysis, we focus on these aspects.
The band structure of nonlinear AM with different nonlinear stiffness (k n ¼ À1 (violet), k n ¼ 0(blue), k n ¼ 1(red)) which are obtained along the M À C À X À M edges of the first Brillouin zone is shown in Fig. 3, where the parameter e¼0:05, a, b, c and A 1 are equal to 0.3. As shown in Fig. 3a, we can observe that due to the existence of the resonator, the band structure of the three stiffnesses will produce a local resonance band gap between the second and third branch. Based on the linear dispersion (blue), as changing nonlinear stiffness, it can be observed that the band structure will also change accordingly. Interestingly, the study found that only the second and third branches have obvious changes within a certain range, which shows that the frequency correction value of the nonlinear stiffness is larger in this area, while the correction value of other areas is smaller. As for the k n ¼ À1 (violet), the frequencies of the first and second branches reduce relatively to linear case, especially in the C À X À M area. It means that the lower boundary of local resonance band gap decreases at k n ¼ À1. Figure 3c shows the details of the third and fourth branches. The frequency of the third branch also has a reduction around the point C, which represents the upper boundary of local resonance band gap. In the case of k n ¼ 1 (red), the opposite phenomenon of k n ¼ À1 occurs, the changing area of the two cases is basically the same, the trend of change in this case is that the frequency of the second-and third-order branches increases. We can find that in the second branch, the point X becomes the lower boundary of band gap, while the upper boundary in the point C has no obvious change in the third branch. It will narrow the band gap. Furthermore, the distributions of band gap for changing nonlinear stiffness k n from À1 to 1 are investigated, as shown in Fig. 3b. The red line represents the lower limit of band gap, the blue line is the upper limit of band gap, and the black line denotes the width of band gap, respectively. It was indicated that with changing nonlinear stiffness from À1 to 1, the frequency of the lower boundary of the band gap gradually rises, while the frequency of the upper boundary increases slightly, resulting in reduction of the width of the band gap. This phenomenon shows that the dispersion relationship is influenced by the nonlinear stiffness k n , and the nonlinear band gap of nonlinear 2D AM can be controlled by adjusting the nonlinear stiffness.

Analysis of wave amplitude range
Let us consider the details of Eq. (30) illustrate the nonlinear dispersion relationship. It demonstrates the importance of the nonlinear resonator to determine the value of X 1 . When the value of X 2 0 ¼ x 2 x 2 n gradually approaches to bc¼ k 3 =m 3 x 2 n , which indicates that x approaches the resonance frequency x 2 re ¼ k 3 m 3 , the value of X 1 would increase and become larger than the value of X 0 and tend to infinity which does not meet the assumptions of the perturbation method used here.
The perturbation analysis assumes an asymptotic development of the term X 1 to obtain the correct result of frequency. According to Eq. (22), we can get the asymptotic development of term To ensure Eq. (27) is true, the asymptotic development of term X 2 should be satisfied as 2eX 0 X 1 ( X 2 0 , and submit the expression of X 1 , we can get the relationship: Because the assuming perturbation parameter e is small, the order of magnitude of ( is equal to e. Then we can get the relationship between wave amplitude A 1 and nonlinear stiffness k n : According to Eq. (32), for a given k n , we can get a set of values of 1 C k n ð Þ in k-space, which represent the allowed wave amplitude for the corresponding wave vector k x ; k y À Á . The result of k n ¼ 1 is plotted in Fig. 4. For convenience, we select the M À C area for analysis. Taking the result of the 1st branch (black) as an example to illustrate this figure, the allowed wave amplitude of the first branch increases as wave number changes from M point to C point and reaches the minimum value around 0.45 at M point. The result of the 2nd branch (red) shows the similar trend, and the minimum value of allowed wave amplitude is around 0.3. Particularly, the value of allowed wave amplitudes of these two branches becomes infinite after the wave number reaches to A and B points, respectively. This phenomenon indicates that almost any value of wave amplitude would not affect the dispersion in the nonlinear regime in the areas A À C and B À C. And the linearization assumption in these two areas remains true even for very high values of amplitude, respectively. This phenomenon can also be observed in Fig. 3a, it is obvious to see that the 1st and 2nd branches change by different nonlinear stiffness k n in the M À A and M À B areas, while the branches have almost no change compared with linear band structure in the areas A À C and B À C. However, as for the cases of the 3rd (blue) and 4th (green) branches, the 3rd branch gets the maximum value 1.75 at point M, while gets the minimum value around 0.55 at point C.
And we can observe that the 4th branch has the same minimum value 0.55 and its maximum value 2.45 was got between the points M and A. In order to analyze the influence of nonlinear effect on the four branches comprehensively, we determine the range of wave amplitude as 0-0.3 (shadow area), because the 0.3 is the maximum of weakly nonlinear system under small parameter assumption. And the nonlinear stiffness k n ¼ 1 is chosen to investigate the amplitude-dependent characteristics of dispersion relationship in the following analysis unless otherwise stated.

Amplitude-dependent dispersion relationship
From the above study, we have obtained the calculation formula of the nonlinear dispersion relationship. It is observed that there is an important parameter wave amplitude A 1 which is quadratic in Eq. (30), and it is necessary to discuss the influence of wave amplitude A 1 on the nonlinear band structures. Figure 5a shows . It can be found that the local resonance band gap exists between the second and third branch and the width would reduce as the wave amplitude A 1 increases. We can also find that the effect of wave amplitude that mostly occurs in local areas (around point X and M À C) is the same as the trend of nonlinear stiffness in Fig. 3. Interestingly, there is a frequency range that the second branch has both with a positive slope and a negative slope around the point X coursed by increasing wave amplitude. This frequency range, called the ''dual wavevector'' region [49], is found in one-dimensional nonlinear elastic metamaterial. Note that at this frequency range, there can be a nonlinear interaction between two waves which have different wavevectors, which causes a complicated wave motion. Figure 5c, d presents the details about the obvious effect of wave amplitude on the 1st and 2nd branches in M À C and C À M areas. We can find that the correction frequency X 1 would increase as the wave amplitude increases, especially around point X. It directly affects the width of the band gap. Figure 5b shows the distributions of the band gap at varying wave amplitude A 1 from 0 to 0.3. The maximum frequency of the second branch (red) gradually increases as wave amplitude increases, while the minimum frequency of the third branch just slightly increases with wave amplitude increases. As a result, the width of band gap shows downward trend as increasing wave amplitude. It can be observed that the wave amplitude plays an important role in nonlinear dispersion relationship of nonlinear 2D AM.

Amplitude-dependent iso-frequency contours and group velocity
The obtained dispersion curves can be used to evaluate the wave properties of the considered nonlinear 2D AM. The group velocity provides important information regarding the directional energy flow and the anisotropic nature of the periodic structures at a specified frequency. The dependence of dispersion on wave amplitude directly affects the group velocity, which suggests how the occurrence of frequencies and directions of preferential propagation may be affected by the amplitude of wave motion. To further illustrate the wave propagation in nonlinear 2D AM as controlled by the wave amplitude, we study the group velocities of the considered nonlinear system at several specific normalized frequencies, on a set of values of wave amplitude A 1 .
According to the definition of group velocity, the group velocity in linear system can be written as: where e Ã denotes reciprocal lattice vectors. Within the perturbation approach described in this study, the group velocity in nonlinear system corresponding to the j À th dispersion can be expressed as: The group velocity is evaluated by first determining the contour of the selected frequency on the dispersion surface, and the gradient of this contour provides the components of the group velocity. From the above research, we discover that the wave amplitude mainly affects the second branch and we will focus on it. Figure 6 shows the three-dimensional representation of iso-frequency contours of the second branch at different wave amplitudes for linear case A 1 ¼ 0 and nonlinear cases A 1 ¼ 0:1; 0:2; 0:3,e¼0:05. From an overall point of view, the three-dimensional representation of iso-frequency contours in the four cases all present a shape of cone. In detail, we can find that all the cases maintain a nearly circular cone shape in the low-frequency range, which shows that linear and nonlinear structures exhibit isotropy in low-frequency range. In this case, the energy flow and directions of wave propagation are mainly to all directions. As for high-frequency range, Fig. 6a shows the linear case. A two-dimensional frequency curve plot at the bottom shows that the structures present strong anisotropy in high-frequency range. With the addition of nonlinearity, it can be found that the wave amplitude mainly affects the iso-frequency contours in high-frequency range. From Fig. 6b-d, we can discover that the iso- frequency contour of the highest frequency will rise with the increase of wave amplitude and gradually form peaks. It also shows strong anisotropy in highfrequency range. This phenomenon implies the amplitude dependent of iso-frequency contours in the nonlinear 2D AM.
Studying wave propagation in periodic structures by using the Bloch theorem with restricting the analysis to a single unit cell is known to all. Usually, it is sufficient to inspect wave vectors lying in the first Brillouin zone. However, in many cases, we only focus on the band gaps, which depend on the extreme values of band structure. Previous research findings are that the extreme value is likely to appear at the boundary of the first Brillouin zone for most symmetric structures. Some works [50,51] have been done to investigate the location of the extreme values of band structure. As for this study, we show the result of dispersion in the whole Brillouin zone to describe the location of extreme values. Figure 7a, b is the twodimensional representation of iso-frequency contour maps of the second branch and the third branch. We plot the profiles of three-dimensional iso-frequency contour maps (as shown in Fig. 6a) corresponding to the blue and red lines in Fig. 7c. The lines denote the central planes parallel to k x À X plane. It can be found that the extreme values of dispersion exist in the contour of Brillouin zone. And it is sufficient to inspect wave vectors in contour of the first Brillouin zone for our structure.
The results for the group velocities corresponding to the second branch for different wave amplitude .') at desired frequency X ¼ 0:4,X ¼ 0:875,X ¼ 0:92 are shown in Fig. 8. As shown in Fig. 8a, at low frequency X ¼ 0:4, we can find that the group velocities of four cases overlap with each other which means that the wave amplitude has almost no effect on group velocity. We can find that the shapes of group velocities present as square and their maximum values are located at diagonal direction that denote the direction of wave propagation. The result also can be clarified in Fig. 6, and the bottom of figures denotes the contour of iso-frequency. As well known, the direction of the external normal of isofrequency curve is the wave propagation direction. Obviously, we can clearly observe its external normal almost along the directions of 45°, 135°, 225°, 315°at the lower frequency. Increasing frequency to X ¼ 0:875, as shown in Fig. 8b, we can find that the shapes of group velocities have a change. And the maximum values of group velocities tend to along the directions of 0°,90°,180°,270°. The wave amplitude affects the amplitude of group velocity that the value of group velocity increases with increasing wave amplitude. As for X ¼ 0:92, shown in Fig. 8c, we can find that the group velocities shape as a cross at small wave amplitude A 1 ¼ 0 and A 1 ¼ 0:1 which shows strong anisotropy, while the group velocities of the larger wave amplitude A 1 ¼ 0:2 and A 1 ¼ 0:3 show a more rounded shape. Compared with Figs. 8a-c, we can find that as the frequency increases, the anisotropy of the group velocity becomes stronger. However, the addition of nonlinearity and changing of wave amplitude will restrain the anisotropy at the same frequency. This is a performance of amplitude-dependent group velocity.

Numerical simulation validation for amplitudedependent dispersion
In this section, we use the time-domain finite-element simulations method [52] to obtain the apparent wave number of the structure to confirm the accuracy of analytical perturbation. The basic principle of this method is to release the system by specifying the initial displacement under corresponding wave (a) (b) (c) Fig. 8 Group velocities of the nonlinear 2D AM for different wave amplitude A 1 at selected frequency X: a number, then observe its evolution with a certain time and finally obtain the frequency of the corresponding wave number. It is noticed that the numerical simulation method is only applicable in theory, because it is unrealistic to directly give a certain initial shape of structure. In this part, we mainly focus on the change of band gap affected by wave amplitude and the value of e. As for dual wavevector, it is hard to observe through numerical simulation [49]. The finite structure represented by an assembly of 60 9 60 unit cells is selected as the research object with unit parameters The results of the linear system under free boundary conditions need to be solved firstly to obtain the values of X 0 from linear dispersion. And the wave amplitude A 3 of the resonator is given by the relationship: It can be used to set the initial condition of the displacement for resonators. Next, we can add the initial conditions to the whole structure: where U rðn;pÞ and V rðn;pÞ denote the principal and resonating displacement of the unit cell located in ðn; pÞ, respectively. The unit cell of 30; 30 ð Þis selected to observe the output signal. In general, in order to avoid the reflection of residual wave, the detection unit cell should be located at the central part of the structure as far as possible. Equation (36) indicates that the initial displacement of mass has a phase shift compared with the unit cell 30; 30 ð Þwithout affecting the results.
Then, given an impose wave number k and collecting the time-domain signal from 30; 30 ð Þ unit cell, the time is selected as t ¼ 20T 0 , where T 0 ¼ 2p X 0 ðkÞ . According to Eq. (36), for example, Fig. 9 shows the initial principal displacements of the unit cells for different wave number k ¼ ½0:5; 0 and k ¼ ½2:72; 0, where n = [1, 60], p = 30. We can find that sinusoidal harmonics appears at k ¼ ½0:5; 0, while space aliasing occurs when wave number increases to 2.72. According to the Shannon theorem, the necessary condition to avoid space aliasing is at least two units per wavelength. In general, sinusoidal harmonics are usually required at six elements per wavelength, and that condition cannot be reached when the wave number exceeds 1. Figure 10 shows one example of output The simulation result of nonlinear dispersion for A 1 is plotted in Fig. 11. Figure 11a shows the nonlinear dispersion for A 1 ¼0:3, e¼0:05, and it is observed that the numerical and theoretical results provide a close match. For convenience, Fig. 11b gives the window enlargement for the first two branches in M À C area. In detail, the results close to the area of C point have a better fit, which shows that nonlinear effects in this part are small. Figure 11c describes the nonlinear dispersion for A 1 ¼0:5 and corresponding simulation result. There is an obvious deviation between the results. According to Sect. 2.2.2, the wave amplitude has a certain range for every branch under the assumption of weak nonlinearity. It is obvious to see that if the given amplitude is larger than the allowed wave amplitude corresponding to the wave vector, the theoretical frequency shift will be larger than the simulation results.
One key point of the weak nonlinear assumption in the perturbation method is the effective range of e, and Fig. 12a shows the frequency shift eX 1 as a function of the e, where A 1 ¼ 0:3, k ¼ ½2:72; 0. We can find a good agreement between theoretical and simulations frequency shift in the range of e ¼ ð0; 0:09Þ. However, simulations results gradually diverge from the theoretical results when e [ 0:09. In detail, the shift frequency results between the theory and simulation have a change of 3:9% at e ¼ 0:09, while the value reaches to 5:4% at e ¼ 0:1. Here, we use the value of 5% as the threshold to measure the failure of results. The value e ¼ 0:09 is almost the upper limit for weak nonlinearity in nonlinear 2D AM, which is close to the previous work [53]. We also note that the theoretical frequency shift is larger than the simulations results at the same e. Wave amplitude is also a critical parameter that affects the weak nonlinearity hypothesis. Figure 12b gives the verification for the effective range of wave amplitude which is proposed in Fig. 4. It is observed that frequency shift of the theoretical and simulations is consistent before around A 1 ¼ 0:3. With increasing wave amplitude, the two calculation results are inconsistent. The nonlinearity of system gradually increases as increasing e and A 1 , and the above research results show that frequency shift will be less than perturbation analysis results when weakly nonlinear assumptions are no longer valid.

Vibration response of finite nonlinear 2D AM
In this section, we will investigate the vibration characteristics of the nonlinear system. The harmonic balance method is introduced to calculate the vibration response of finite nonlinear 2D AM.
Let us consider the finite nonlinear 2D AM with 5 9 5 unit cells as shown in Fig. 1.
For the wave propagation along the x-direction, the equations of motion are given by where U ¼ ½u 1 1;1 ð Þ q 1 1;1 ð Þ u 1 1;2 ð Þ q 1 1;2 ð Þ . . .u 1 n;p ð Þ q 1 n;p ð Þ T is the displacement vector, which contains the displacement q 1 n;p ð Þ ¼ u 3ðn;pÞ À u 1ðn;pÞ . M denotes the mass matrix of finite 2D AM, and K and K N ðk n Þ are the linear and nonlinear stiffness matrix of the nonlinear system. F ¼ ½f cosðxtÞ 0 . . . 0 T is the force vector. There is only one harmonic excitation force acting on the (1,1)th unit cell.
Applying the harmonic balance method to solve frequency response, the solution of the system is assumed as: Submitting Eq. (38-37), we can obtain: then using the trigonometry relationships (sin 3 xt ð Þ ¼ 3 4 sin xt ð Þ þ 1 4 sin 3xt ð Þ,-) and ignoring the high-order harmonic terms, we can get: Equation (40) can be solved and the response amplitude under the excitation force is Y ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi a 2 þ b 2 p . Finally, we can get the frequency response T is defined by: where the Y out and Y inp are the amplitudes of output and excitation point, respectively. Usually, we select the unit cells of left boundary that the (1,1)th to (5,1)th unit cell is the excitation point, while unit cells of right boundary are the output point. It should be noted that due to the different methods for nonlinear system, the parameters in this section are not exactly equal to Sect. 2, but the trend of vibration characteristics affected by nonlinear effects should follow the same way. Figure 13 shows the frequency response with different excitation amplitude f , where f ¼ 5ðredÞ, f ¼ 10ðblueÞ, k n ¼ 1. The comparison of the results in Fig. 13a, b shows that increasing the number of the cells will improve the accuracy of results. We can find that the frequency response T of a linear system (black) is independent of the excitation amplitude, while the excitation amplitude has an important influence on the nonlinear system. It is obvious to see that a strong attenuation exists between X ¼ 1 and X ¼ 1:5 which present the band gap that elastic waves cannot propagate. As the excitation amplitude f is increasing, the frequency of lower boundary of attenuation area will increase slightly. This excitation amplitude-dependent behavior of nonlinear 2D AM is similar to the trend of wave amplitude-dependent band structure in the previous study (Fig. 5).

Conclusions
In this study, we consider a type of nonlinear twodimensional acoustic metamaterials. The nonlinear systems are made of basic linear two-dimensional spring-mass system and distributed Duffing resonators. The perturbation method is introduced to evaluate the nonlinear dispersion relationships and wave propagation characteristics. A quantitative study of the maximum admissible amplitude has been performed to provide an overview of the limits of the proposed method. Furthermore, the harmonic balance method is used to investigate the frequency response and effective mass of the nonlinear systems.
The results show that the two-dimensional nonlinear acoustic metamaterials can tune the characteristics of wave propagation and effective mass in the structures by adjusting the nonlinear stiffness and wave amplitude within a certain range. The main findings are listed below: (1) Nonlinear stiffness and wave amplitude can significantly affect the dispersion curves. The maximum allowable wave amplitude was found under the assumption of small parameters. The influence of nonlinear stiffness and wave amplitude on the dispersion curve is similar. They mainly affect the second branch of band structure, while have a little effect on other branches. The distributions of band gap indicate that increasing the nonlinear stiffness and wave amplitude will narrow the band gap. Numerical simulations validate the dispersion results and ensure the effective range of the value of e. (2) The iso-frequency contours and group velocity of nonlinear 2D AM both present amplitudedependent. The iso-frequency contour of the highest frequency will rise with the increase of wave amplitude and gradually form peaks which represent the anisotropic behavior, while in the group velocity plots, increasing wave amplitudes will restrain the anisotropic behavior at the same frequency. (3) The vibration response and effective mass of nonlinear 2D AM also present the amplitudedependent. And the area of wave attenuation moves to higher frequency as increasing excitation amplitude. There is a negative mass region, and its frequency range moves to higher frequency as increasing wave amplitude.

Declarations
Conflict of interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.