## Abstract

Wave propagations exhibit direction and frequency selectivity in two-dimensional (2D) periodic structures, which provides possibilities to regulate wave dispersion and bandgap properties. Most of current researches focus on regulations of 1D waves, and there are few works about active regulations of 2D waves, especially in the structures with strong nonlinearities that have remarkable influences on dispersions. In this work, two types of 2D periodic nonlinear lattice structures with piezoelectric springs, which include a monatomic and a diatomic structure, are designed to implement controllable dispersion and propagation direction of 2D waves. Considering the strong nonlinearities caused by the cubic spring, dynamic models of the wave propagations in the two kinds of periodic structures are established, and an improved incremental harmonic balance (IHB) method is developed to implement efficient and accurate calculations of the 2D wave propagation. Influences of active and structural parameters on dispersion and bandgap properties are comprehensively studied, and the regulation ability of the piezoelectric springs is demonstrated where the proportional voltage constant is the active control parameter with particle displacements as the feedback. Results also show that a piezoelectric modulated bandgap and a critical wave vector region are created by positive and negative proportional constants, respectively, which indicate that the structures can be used to filter a wide range of low-frequency long-wavelength noises and waves at particular directions. The properties predicted by the improved IHB method are verified by numerical experiments.

## 1 Introduction

Periodic materials and structures have been extensively studied in wave propagation problems over the past decades [1]. The periodic structures can exhibit bandgap properties that prohibit wave propagations in one or more directions [2]. The bandgap properties of periodic structures also have promising applications in vibration absorption, noise reduction, and acoustic modulation [3–5]. To improve the bandgap characteristics, researchers have conducted extensive studies on different periodic structures [6–8].

Adjustable bandgap properties usually include the number, position, and width of the bandgaps [9–12]. The bandgap characteristics of periodic structures that work in linear regions are designed by inclusions, geometry, and elastic characteristics [13]. Structures that work at nonlinear regions can enhance bandgap sensitivity and improve the adjustability of the structures [14–16]. At present, 1D nonlinear periodic structures are widely used as theoretical models in wave propagation problems [15,17,18]. In 2D periodic structures, not only dispersions but also directionality characteristics are shown in wave propagations, which provides new methods to implement wave beaming, spatial filtering, and imaging [19,20].

In the above-mentioned studies, bandgap characteristics are passively controlled by structure properties, which are unchangeable in use. Active regulation of the bandgap characteristics can be adaptive to environments and have a larger range of adjustment [21,22]. There are two mechanisms of active regulation. One is mechanical reconfigurability [23]. The other mechanism is to apply multi-physics fields [24–26]. Piezoelectric materials are widely used in current research [27]. And piezoelectric springs have been successfully used for bandgap regulations of 1D periodic structures [28,29]. However, most of the mentioned research studies mainly focus on regulations of 1D waves, and there are few research studies on active control of 2D nonlinear structures.

Furthermore, the mechanisms of bandgap adjustment have to be comprehensively studied by exploring the influences of structural and active electric parameters on bandgaps. There is a need to develop an effective computational method to obtain dispersion relationships for structures with general nonlinearities. Currently, the perturbation method is often used to solve 2D-lattice wave problems with weak nonlinearities [30]. The incremental harmonic balance (IHB) method can be used to solve wave problems with strong nonlinearities. To avoid complex and time-consuming calculations, low-order trigonometric functions are usually utilized as the basis functions [30], which will decrease the accuracy when a strong nonlinearity is involved. The IHB method was improved in studies of 1D strong nonlinear wave problems [31–34]. By understanding the similarity between steady-state wave and vibration, the 1D wave problem can be transformed into a delay differential equation (DDE), which can avoid redundant basis functions in the description of wave responses, and both the accuracy and efficiency can be improved [33]. However, the current improved IHB method for 1D structures cannot be applied directly to the 2D problems. Thus, an algorithm framework of an improved IHB method is highly required to deal with 2 or high-dimensional wave propagation problems.

This work focuses on active bandgap regulations on 2D periodic structures, and an efficient semi-analytical method is derived to calculate 2D wave propagation with general nonlinearities. The article consists of six parts. In the next section, a type of 2D nonlinear tunable periodic structure is constructed, and dynamic equations are derived. In Sec. 3, an improved IHB method for 2D wave problems is derived, and its accuracy and efficiency are verified. In Sec. 4, dispersion and bandgap properties of 2D waves in the structures are comprehensively studied with active and structural parameters for regulations. In Sec. 5, the dispersion properties obtained by the improved IHB method are verified by numerical experiments. Finally, results of the work are summarized.

## 2 2D Periodic Structures With Piezoelectric Springs

### 2.1 Definition of a 2D Periodic Structure.

2D Polyatomic lattice systems with periodic arrangement are studied, and a general arrangement model is represented in the *X*–*Y* plane, as shown in Fig. 1. Because of periodicity, the arrangement of each type of atom is the same, and there is a global translation between different types of atoms. Vectors *a*_{1} and *a*_{2} are formed by two adjacent atoms that belong to the same type, and the atoms are arranged periodically along the vectors *a*_{1} and *a*_{2}. In the quadrilateral formed by vectors *a*_{1} and *a*_{2}, vertexes represent lumped masses of atoms, and edges represent interactions between adjacent atoms, which have nonlinear interactions in this work. Integers *n*_{1} and *n*_{2} represent locations of atoms, and origin O represents the point where *n*_{1} and *n*_{2} are both zeros.

**represent the wave vector in the inverted lattice space**

*k**μ*

_{1}and

*μ*

_{2}represent wave numbers in directions of

*b*

_{1}and

*b*

_{2},

*b*

_{1}and

*b*

_{2}represent the corresponding vectors of

*a*

_{1}and

*a*

_{2}in the inverted lattice space, which are perpendicular to

*a*

_{1}and

*a*

_{2}, respectively.

*a*)

*b*)

*c*)

### 2.2 Tunable 2D Monatomic Lattice.

An actively tunable 2D monatomic lattice structure is designed, which consists of atoms of an equal mass *m* arranged in the form of a square lattice on the *X*–*Y* plane, as shown in Fig. 2(a). Each atom is connected by five springs, which are four in-plane springs with stiffness *k _{x}* and

*k*in

_{y}*X*- and

*Y*-directions, respectively, and a piezoelectric spring with stiffness

*k*in the

_{E}*Z*-direction, as shown in Fig. 2(b). The other end of the piezoelectric spring is connected to the base, and the control voltage can be applied through the base.

Transverse wave propagation in the *X*–*Y* plane is studied, so that each particle vibrates in the Z-direction. A schematic diagram of the 2D structure model is shown in Fig. 3. Vectors *a*_{1} and *a*_{2} are the propagation directions of the wave with wave numbers *μ*_{1} and *μ*_{2}, respectively (Fig. 3). The first Brillouin region in the inverted lattice space of the structure is shown in Fig. 3, where the high symmetric points are Γ = (0, 0), *X* = (0, *π*), *M* = (*π*, *π*). The subsequent studies on dispersion properties are mainly carried out along the boundary of the region.

Since each atom is constrained to move only in the *Z*-direction, *Z*-components of the five spring forces are developed here. The force of the piezoelectric spring is directly in the *Z*-direction. Forces of the *X*- and *Y*-direction springs are studied in the *X*–*Z* and *Y*–*Z* planes, respectively, as shown in Fig. 4. The initial length of the spring is *a*, *F _{kx}* and

*F*represent forces of the

_{ky}*X*- and

*Y*-direction springs, respectively, along the spring direction, and

*F*and

_{zx}*F*represent

_{zy}*Z*-components of

*F*and

_{kx}*F*, respectively. If there are pretensions in the springs, the initial pre-stretch of the springs in both directions is the same,

_{ky}*δ*. The

*Z*-direction displacement of the atom at (

*n*

_{1},

*n*

_{2}) is $u(n1,n2)$.

*n*

_{1},

*n*

_{2}) and (

*n*

_{1}−1,

*n*

_{2}) can be expressed as

*Z*-direction can be expressed as

*Z*-direction force component of the

*Y*-direction spring is

It can be seen that the in-plane springs exhibits nonlinear characteristics in the *Z*-direction, and its linear stiffness can be controlled by the pretension.

*n*

_{1},

*n*

_{2}) can be expressed as

*k*represents the stiffness of the piezoelectric spring and its equilibrium is zero. Based on the constitutive equations of the piezoelectric spring (Appendix A), one has

_{E}*h*is the piezoelectric coupling constant,

_{p}*κ*is the dielectric constant,

^{S}*C*is the elastic modulus,

^{D}*b*,

*t*and

_{p}*L*are the width, thickness, and length of the piezoelectric spring, respectively, and

_{p}*V*is the applied voltage, which can be controlled by the external circuit. In the work, a linear control law based on displacement feedback is designed to control piezoelectric voltage

_{p}*K*represents the proportional control constant. Combining Eqs. (11) and (12), the stiffness of piezoelectric spring in Eq. (10) can be expressed as

_{g}*a*)

*b*)

*c*)

*k*and

_{e}*k*are the active and structural stiffness of the piezoelectric spring, respectively. It can be seen that the stiffness of the piezoelectric spring can be regulated by the proportional constant. By combining Eqs. (9) and (10), the dynamic equation of the atom at (

_{s}*n*

_{1},

*n*

_{2}) is

### 2.3 Tunable 2D Diatomic Lattice.

The tunable 2D monatomic structure presented in Sec. 2.2 can be further extended to multi-atomic structures. In this section, an actively tunable 2D diatomic lattice structure is designed, which consists of atoms of different masses *m*_{1} and *m*_{2} in the form of a square lattice on the *X*–*Y* plane, as shown in Fig. 5(a). Each atom is connected by five springs, which are four in-plane springs with stiffness *k*_{x} and *k*_{y} in X- and Y-directions, respectively, and a piezoelectric spring with stiffness *k _{E}* in the

*Z*-direction, as shown in Fig. 5(b). The other end of the piezoelectric spring is connected to the base, and the control voltage can be applied through the base. The displacement of the atoms with masses

*m*

_{1}and

*m*

_{2}can be expressed as

*u*

_{1}and

*u*

_{2}, respectively.

Transverse wave propagations of the diatomic lattice are studied in the *X*–*Y* plane, and each particle vibrates in the *Z*-direction. A schematic diagram of the tunable 2D diatomic structure is shown in Fig. 6. Vectors *a*_{1} and *a*_{2} are the propagation directions of the wave with wave numbers *μ*_{1} and *μ*_{2}, respectively. The first Brillouin region in the inverted lattice space of the structure is shown in Fig. 6, where the high symmetric points are $\Gamma =(0,0),X=(\pi 4,\pi 4),M=(0,\pi 2)$. The subsequent studies on dispersion properties are mainly carried out along the boundary of the region.

*n*

_{1},

*n*

_{2}) can be expressed as

*u*

_{1}and

*u*

_{2}at (

*n*

_{1},

*n*

_{2}), respectively, and

*a*)

*b*)

To sum up, both monatomic and diatomic lattice structures constructed in the section are nonlinear model structures, and the intensity of the nonlinearity is related to the parameters in the governing equations. Therefore, an efficient analysis method suitable for general nonlinear systems is required. An improved IHB method for 2D wave problems is introduced as follows.

## 3 Solving 2D Nonlinear Wave Problems by an Improved Incremental Harmonic Balance Method

### 3.1 Method Derivation.

*F*(·) is the nonlinear function describing wave propagations,

*t*is the time variable, and the overdot represents the time derivative, and $U(n1,n2)$ represents particle motions of the unit cell at position (

*n*

_{1},

*n*

_{2}). For the transverse wave, each particle has one-degree-of-freedom and moves in the out-of-plane direction:

*k*th particle in the cell at (

*n*

_{1},

*n*

_{2}), with

*k*= 1, …,

*n*. The function

*F*(·) is

*n*-dimensional function

*F*= [

*F*

_{1},

*F*

_{2},…,

*F*

_{n}]

^{T}, where each component represents the governing equation of a particle in the cell.

*n*

_{1},

*n*

_{2}) = (0, 0), Eq. (17) can be written as

*i*,

*j*) can be expressed as

*ω*represents the angular frequency, and

*T*

_{(i,j)}represents the wave vector of the unit

*μ*

_{1}and

*μ*

_{2}represent wave numbers in the 2D medium.

*U*(

*ωt*−

*T*) =

*C*

_{m}

*Q*, where

*C*

_{m}is a matrix of truncated Fourier basis functions and

*Q*is a matrix of coefficient vectors that are going to be determined. By constructing the first-derivative, second-derivative, and time-delay operators, i.e.,

*G*

_{m}_{1},

*G*

_{m}_{2}, and $Dm(i,j)$, respectively, the dynamic equation in Eq. (22) can be converted to an algebraic equation with the underdetermined variable

*Q*:

The derivation process is shown in Appendix B.

*Q*in Eq. (23):

*R*is the harmonic residual of the function

*F*in Eq. (23) with a guess of the variables

*Q*and

*ω*,

*J*and

_{Q}*J*represent Jacobian matrices of

_{ω}*Q*and

*ω*, and Δ

*Q*and Δ

*ω*represent increment of the variables, respectively. The derivation of the residual is shown in Appendix C, and that of the semi-analytic form of the Jacobian matrix

*J*is shown in Appendix D:

_{Q}The Jacobian matrix *J _{ω}* can be obtained by the Fourier transform of $\u2202F(CmGm2Q,\tau ,\omega )/\u2202\omega $.

### 3.2 Method Verification.

*ε*is used to adjust the strength of the nonlinear term since the nonlinearity is weak in the perturbation method [30]. The residual and Jacobian matrices of Eq. (26) with the use of the improved IHB method are given in Appendix E. Since there is no piezoelectric spring in Ref. [30], the stiffness of the piezoelectric spring in Eq. (26) is set to be zero (i.e.,

*β*

_{s}=

*β*

_{e}= 0), and the system in Fig. 2 degenerates to a conventional cubic-spring model.

The reference parameters are defined as $m=1kg,k1=1Nm\u22121,$$k2=1.5Nm\u22121,\Gamma 1=\Gamma 2=1Nm\u22121$. When $\epsilon =0.01,A=0.002$, the system has a weak nonlinearity, and the results of the improved IHB method are compared with those of the perturbation method, as shown in Fig. 7(a). The dispersion curves obtained by the two methods are exactly the same, which verifies the accuracy of the improved IHB method for 2D wave problems with weak nonlinearities. When $\epsilon =1,A=2$ and the reference parameters are unchanged, the system has a strong nonlinearity, and the dispersion curves calculated from the improved and original IHB methods are shown in Fig. 7(b), where the Fourier expansion order is nine. Results show that the curves match well and the improved IHB method can be used to solve 2D wave problems with strong nonlinearities.

Since the Fourier expansion order of the IHB methods determines the accuracy of the solution, the convergence analysis of the order is needed. The response frequencies of steady-state waves with different expansion orders *M* are shown in Table 1. The frequencies obtained by low-order basis functions are different from those obtained by high-order basis functions for some wave vectors, indicating that the results do not converge at low expansion orders. For example, when ** κ** = (0.5, 0.5), the frequency with the first-order expansion is 0.9876, and those of the seventh- and ninth-order expansions are both 1.4046. The results indicate that the accuracy cannot be guaranteed with the use of the first-order expansion in Ref. [31]. When the expansion order

*M*≥ 7, the convergence of steady-state wave responses is guaranteed in the first Brillouin region, and

*M*= 7 is used in the following calculation.

κ | (0.5, 0.5) | (1, 1) | (1.5, 1.5) | (2, 2) | (2.5, 2.5) |
---|---|---|---|---|---|

M = 1 | 0.9876 | 2.7241 | 5.0575 | 7.4628 | 9.3618 |

M = 3 | 1.4962 | 3.7463 | 5.1712 | 7.4658 | 9.4171 |

M = 5 | 1.3082 | 3.0835 | 5.1720 | 7.4658 | 9.4171 |

M = 7,9 | 1.4046 | 3.0837 | 5.1720 | 7.4658 | 9.4171 |

κ | (0.5, 0.5) | (1, 1) | (1.5, 1.5) | (2, 2) | (2.5, 2.5) |
---|---|---|---|---|---|

M = 1 | 0.9876 | 2.7241 | 5.0575 | 7.4628 | 9.3618 |

M = 3 | 1.4962 | 3.7463 | 5.1712 | 7.4658 | 9.4171 |

M = 5 | 1.3082 | 3.0835 | 5.1720 | 7.4658 | 9.4171 |

M = 7,9 | 1.4046 | 3.0837 | 5.1720 | 7.4658 | 9.4171 |

The calculation efficiency of the method is studied. In the calculation process, the method does not need to carry out repeated derivation work, which greatly simplifies the use process of a new problem. The calculation time of two IHB methods coded in matlab is shown in Table 2 to demonstrate the efficiency. When the wave vector ** κ** = (0.5, 0.5), the time of the improved IHB method is always less than 1

*s*, and that of the IHB method [58] can be very long when expansion orders are high. To study the dispersion property of the system, a large number of wave vectors are selected, and a method with high calculation efficiency is expected.

## 4 Active Regulations of Bandgap Properties

### 4.1 Bandgap Regulations of the 2D Monatomic Structure.

Regulation mechanisms of bandgap properties are studied for the 2D monatomic lattice structure designed in Sec. 2.2, where electric voltages are the active regulation parameters and influences of structural parameters are also studied. The reference parameters in Eq. (14) are set as $m=1kg,k1=k2=1Nm\u22121,\Gamma 1=\Gamma 2=1Nm\u22121,A=0.002$, and dispersion properties with different dimensionless stiffness of the piezoelectric spring *β*_{s} and *β*_{e} are shown in Fig. 8. The active stiffness *β*_{e} can be changed from negative to positive by the proportional constant *K _{g}*, as shown in Figs. 7(a)–7(c).

When the stiffness *β*_{e} = *β*_{s} = 0, the dispersion curve of the system is the same as that of the cubic-spring model. When *β*_{e} > 0, the system is in positive proportional control (positive *K _{g}*), and a new bandgap appears in the low-frequency range, which is called a piezoelectric bandgap, as shown in Fig. 8(b). The existence of the bandgap means that there are no steady-state motions of particles with the frequencies in the bandgap at the wave vector Γ (0, 0). The reason is that particle motions at Γ (0, 0) are pure vibrations and the resonant frequency is nonzero when the piezoelectric spring stiffness is positive. Thus, low-frequency waves can be isolated by the piezoelectric-spring-controlled structures. When

*β*

_{e}< 0 (negative

*K*), the system is in negative proportional control, and a critical wave vector region exits, as shown in Fig. 8(c). The critical region means that only particle vibrations with sufficient phase differences can be stably propagated in the structure. The reason is that the particle vibrations are not stable with the negative stiffness, and a stable motion only exists if there are sufficient recovery forces provided by adjacent particles, which requires a large enough phase difference between the particle motions. Therefore, the critical wave vector region exists near Γ (0, 0), and waves in the region, which have long wavelengths, cannot be propagated in the structure at the negative control constant.

_{g}The influence of the proportional constant on the system is studied. The width of the piezoelectric bandgap increases with the positive *β*_{e}, as shown in Fig. 8(b), which can be actively controlled by an external circuit. The width of the critical wave vector region increases with the negative *β*_{e}, as shown in Fig. 8(c), which indicates that the directional selectivity of the 2D wave propagation can be actively adjusted. Although the dimensionless structural stiffness *β*_{s} cannot be online adjusted, increasing *β*_{s} yields a more sensitive control of the piezoelectric bandgap width, as shown in Fig. 8(d).

The heatmap of the dimensionless frequency with respect to the wave vector (*μ*_{1}, *μ*_{2}) is shown in Fig. 9. The heatmaps of the zero piezoelectric spring stiffness and positive spring stiffness (*β*_{e} = 0.2, *β*_{s} = 0.1) are generated with wave numbers ranging from 0 to *π*. The results in the first Brillouin region are the same as those in Fig. 9, and frequencies close to Γ (0, 0) are not zero, as shown in Fig. 8(b), indicating the existence of the piezoelectric bandgap. When *μ*_{1} = *μ*_{2}, the wave has the widest frequency range. When *μ*_{1} = 0 or *μ*_{2} = 0, the wave propagates with the smallest frequency range.

The relationship between the dispersion and other structural parameters of the system is studied here, where the dimensionless piezoelectric spring stiffness is defined as *β*_{e} = 0.2, *β*_{s} = 0.1. When the linear spring stiffness in the *Y*-direction increases, the frequency curve of the dispersion rises, as shown in Fig. 10(a). It is worth noting that the segment of the frequency curve in the wave vector interval Γ → *X* does not change with the stiffness ratio *ϕ*, because there is no wave vector in the *Y*-direction. The frequency also increases with nonlinear spring stiffness, and the change is unremarkable near Γ (0, 0), as shown in Figs. 10(b) and 10(c). Results show that the frequency increases with the amplitude-induced nonlinearity, as shown in Fig. 10(d), and the above parameters do not affect the piezoelectric bandgap.

### 4.2 Bandgap Regulations of the 2D Diatomic Structure.

*a*)

*b*)

The residual and Jacobian matrices of Eq. (27) with the use of the improved IHB method are given in Appendix F. The reference parameters in Eq. (27) are defined as $m1=1kg,m2=2kg,k1=$$1Nm\u22121,k2=1Nm\u22121,\Gamma 1=\Gamma 2=1Nm\u22121,andA=0.002$, and dispersion properties with different dimensionless stiffness of the piezoelectric spring *β*_{s} and *β*_{e} are shown in Fig. 10. The active stiffness *β*_{e} can be changed from the negative to positive by the proportional constant *K _{g}*, as shown in Figs. 11(a)–11(c). When the stiffness

*β*

_{e}=

*β*

_{s}= 0, the dispersion curve of the system is the same as that of the cubic-spring model, and there is a structural bandgap between the acoustic and optical wave curves. When

*β*

_{e}> 0, a piezoelectric bandgap of the diatomic structure appears in the low-frequency range, as shown in Fig. 11(b). It is worth noting that a critical wave vector region that only appears in acoustic waves is also created in the diatomic structure by a negative proportional (

*β*

_{e}< 0) constant, and the range is larger than that of the monatomic structure so that the low-frequency waves not only with long wavelengths but also at particular directions cannot be propagated in the diatomic structure, as shown in Fig. 11(c).

The influence of the proportional constant on the system is studied. The width of the piezoelectric bandgap increases with the positive *β*_{e}, as shown in Fig. 11(b), which can be actively controlled by the external circuit. The width of the critical wave vector region increases with the negative *β*_{e}, which can cover part of the sector *X* → *M* and Γ → *M*, as shown in Fig. 11(c). It indicates that the directional selectivity of 2D acoustic waves can be actively adjusted. Although the dimensionless structural stiffness *β*_{s} cannot be online adjusted, increasing *β*_{s} yields a more sensitive control of the piezoelectric bandgap width, as shown in Fig. 11(d), which is similar to that of the monatomic structure. Remarkably, the stiffness of the piezoelectric spring has little effect on the width of the structural bandgap but has obvious effects on the frequency ranges of the dispersion curves.

The relationship between the dispersion and other structural parameters of the system is studied here. In the analysis, the dimensionless piezoelectric spring stiffness is defined as *β*_{e} = 0.2, *β*_{s} = 0.1. When the mass ratio of two particles increases, the width of the piezoelectric bandgap decreases, the width of the structural bandgap increases, and the range of the dispersion reduces, as shown in Fig. 12(a). The reason is that when the wave vector is at Γ (0, 0), the governing equation of the system can be regarded as the vibration equation, and the mass ratio of the particles can regulate the resonance frequency of the system. When the linear spring stiffness in the *Y*-direction increases, widths of the structural and piezoelectric bandgaps are almost unchanged, and the frequency curve of the dispersion rises, as shown in Fig. 12(b). The frequency also increases with nonlinear spring stiffness, and the change is unremarkable near Γ (0, 0), as shown in Fig. 12(c). The frequency increases with the amplitude-induced nonlinearity, as shown in Fig. 11(d), and the parameters do not affect the piezoelectric bandgap.

## 5 Numerical Experiments

To further verify the bandgap properties of the above model, wave propagation results are numerically studied by a simulation model that is developed in matlab, as shown in Fig. 13. In the simulation, an input excitation is exerted on the unit cell at (25, 1), and the motion of the cell at (25, 25) is the output. The fourth-order Runge–Kutta numerical integration method is used to solve the wave problem given in Eqs. (14) and (15). The periodic external excitation is $F=fsin(\omega t)$, and the transient and steady-state responses of the output are recorded to study the wave propagation properties.

The bandgap frequency of the 2D monatomic lattice is investigated, where the system parameters are $m=1kg,ke=0.2Nm\u22121,$$ks=0.1Nm\u22121,k1=1Nm\u22121,k2=1Nm\u22121,\Gamma 1=\Gamma 2=1Nm\u22121$. From the steady-state results obtained by the IHB method (Fig. 8), there is no elastic wave when *ω*_{n} = 0.2 in the piezoelectric bandgap. The numerical simulation shows that the vibration of the particle vanishes at the frequency, as shown in Fig. 14(a), which verifies the predicted phenomenon by the IHB method. When the dimensionless frequency *ω*_{n} = 1, there is an elastic wave by the IHB method, which is also verified by the simulation in Fig. 14(a).

Next, the bandgap frequency of the 2D diatomic lattice is investigated, where the system parameters are $m1=1kg,m2=2kg,$$ke=0.2Nm\u22121,ks=0.1Nm\u22121,k1=1Nm\u22121,k2=1Nm\u22121,\Gamma 1=\Gamma 2=1$$Nm\u22121$. From the steady-state results obtained by the IHB method (Fig. 11), there are no elastic waves when *ω*_{n} = 0.2 and *ω*_{n} = 1.5 in the piezoelectric bandgap. Numerical simulations also show that vibrations of the particle vanish at the frequencies, as shown in Figs. 14(b) and 14(c), which verify the predicted phenomenon by the IHB method. In Fig. 14(c), the particle may have resonated when *ω*_{n} = 1.5, causing the particle to be strongly excited, but the motion of the particle is suppressed due to the action of the bandgap. When *ω*_{n} = 2.2, after the particle is excited, the motion displacement decreases and then increases and tends to be stable, which may be related to the refraction and reflection of the elastic wave when propagating, causing the coherence or cancellation of the elastic wave, resulting in the amplitude of the stability is less than the amplitude of the excitation.

## 6 Conclusion

In this work, a type of 2D nonlinear periodic lattice structures, i.e., the monatomic and diatomic structures, with the control of piezoelectric springs, is designed. With the proportional constant as the control variable, the bandgap properties of the structures can be actively regulated, and the filtering characteristics in specific frequencies and directions can be realized. The proposed design has a simple structure, and it is easy to implement and control. Wave vectors in 2D structures are introduced into the improved IHB method to realize the DDE transformation of the governing equations, which extends the improved IHB method to multi-dimensional problems. Accuracy of the improved IHB method is verified, and computational efficiency is higher than other methods. With the method, active regulations of dispersion and bandgap properties of the established structures are studied with respect to the structural and active parameters of the system. The regulation properties predicted by the IHB method are verified by numerical experiments. Highlights of the regulation properties are as follows:

new bandgaps (piezoelectric bandgaps) of waves in the structures are created by a positive proportional constant, which means that low-frequency waves can be isolated by the piezoelectric-spring-controlled structures;

a critical wave vector region is created in the monatomic structure when the constant is negative, which means that waves with long wavelengths cannot be propagated in the structure with a negative control constant; and

a wider critical wave vector region that only appears in acoustic waves in the diatomic structure is created by the negative proportional constant, which means that low-frequency waves not only with long wavelengths but also at particular directions cannot be propagated in the diatomic structure.

Based on the properties, the structures with the feedback control of the piezoelectric spring can be used to filter a range of low-frequency and long-wavelength noises, as well as those at particular directions.

## Acknowledgment

We gratefully acknowledge support from the National Key R&D Program of China (Grant No. 2022YFB3203600) and the National Natural Science Foundation of China (Grant No. U1930206).

## Funding Data

The National Key R&D Program of China (Grant No. 2022YFB3203600).

The National Natural Science Foundation of China (Grant No. U1930206).

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Appendix A: Constitutive Relation of Piezoelectric Spring

*E*is the electric field intensity,

_{p}*D*is the electric displacement,

_{p}*T*is the stress,

_{p}*S*is the strain,

_{p}*h*is the piezoelectric coupling constant,

_{p}*κ*is the dielectric constant, and

^{S}*C*is the elastic modulus. By introducing some new physical quantities, Eqs. (A1) and (A2) can be expressed as

^{D}*b*,

*t*, and

_{p}*L*are the width, thickness, and length of the piezoelectric spring, respectively;

_{p}*Vp*is the applied voltage,

*F*is the compression force, and

_{p}*Q*is the charge. According to Eqs. (A3) and (A4), it can be expressed as

_{p}## Appendix B: Equivalent Transformation of Equation

*T*represents the transpose of the vector, and

*M*represents the order of Fourier truncation, $Q=[q1T,q2T,\cdots ,qnT]T,qk=[a0k,a1k,\cdots ,aMk,b1k,\cdots ,bMk]T,k=1,\cdots ,n$ is the coefficient vector of

*C*,

_{S}*G*

_{m1}, and

*G*

_{m2}are square matrices of

*n*(2

*M*+ 1)-order

*H*

_{1}and

*H*

_{2}are first and second-derivative matrix operators, respectively:

## Appendix C: Residual Calculation

*R*can be expressed as

*F*:

*C*

_{m}_{2}can be expressed as

*C*

_{s}_{1}is an orthogonal basis for

*C*

_{s}## Appendix D: Semi-Analytic Expression of the Jacobian Matrix

*J*for the 2D wave problem is constructed as

_{Q}*f*

_{DD}can be expressed as

*k*= 1, …,

*n*is the Toeplitz matrix form of the coefficient vector of

*M*th-order Fourier functions of $\u2202Fk\u2202U\xa8$, $[f\u2212Mk\cdots f0k\cdots fMk]$ that can be calculated by fast Fourier transform (FFT):

*f*

_{D}and

*f*

_{(i,j)}can be expressed respectively as

*k*= 1, …,

*n*, and $f(i,j)k1,k2$,

*k*

_{1, 2}= 1, …,

*n*, represent the Toeplitz matrix form of the coefficient vectors of $\u2202Fk\u2202U\u02d9$ and $\u2202Fk\u2202U(i,j)$, respectively.

## Appendix E: Derivation of Monatomic Lattice by Improved IHB Method

## Appendix F: Derivation of Diatomic Lattice by Improved IHB Method

*T*

_{2}= −1 − 3

*λ*

_{1}(

*u*

_{1}−

*u*

_{2})

^{2},

*T*

_{3}= −1 − 3

*λ*

_{1}(

*u*

_{2}−

*u*

_{1})

^{2}

_{,}and $T4=(\beta s+\beta e)+2(1+\varphi )+3\lambda 1((u2\u2212u1)2+(u2\u2212u1(1,0))2)+3\lambda 2((u2\u2212u1(0,1))2+(u2\u2212u1(1,\u22121))2)$