Abstract
Numerous recent works have established the potential of various types of metamaterials for simultaneous vibration control and energy harvesting. In this paper, we investigate a weakly nonlinear metamaterial with electromechanical (EM) local resonators coupled to a resistance-inductance shunt circuit, a system with no previous examination in the literature. An analytical solution is developed for the system, using the perturbation method of multiple scales, and validated through direct numerical integration. The resulting linear and nonlinear band structures are used for parametric analysis of the system, focusing on the effect of resonator and shunt circuit parameters on band gap formation and vibration attenuation. This band structure analysis informs further study of the system through wavepacket excitation as well as spectro-spatial analysis. The voltage response of the system is studied through spatial profiles and spectrograms to observe the effects of shunt inductance, nonlinearity, and their interactions. Results describe the impact of adding a shunted inductor, including significant changes to the band structure; multiple methods of tuning band gaps and pass bands of the system; and changes to wave propagation and voltage response. The results demonstrate the flexibility of the proposed metamaterial and its potential for both vibration control and energy harvesting, specifically compared to a previously studied system with resistance-only shunt.
1 Introduction
Metamaterials are artificially engineered structures that possess properties not found in naturally occurring materials [1]. These properties range from zero or negative values of standard engineering parameters, such as density and Poisson's ratio [2], to nonlinear phenomena, such as gap solitons [3] and asymmetric wave propagation [4]. Metamaterials have a foundation in optics and electromechanics, and include mechanical or acoustic metastructures that exploit elastic and wave properties such as motion, deformations, stresses, and mechanical energy [5–7]. The unusual features of metamaterials make them beneficial for numerous applications, including vibration and noise control, energy harvesting, nondestructive testing, and acoustic rectifiers.
Metamaterials consist of many unit cells arranged in periodic or aperiodic patterns. A periodic design can allow them to display properties and functionalities that exceed those of their components. As a result, the study of periodic structures forms a major part of the foundation for metamaterial research [8–12]. In these and other works, it was observed that periodic structures prevent waves from propagating through the structure at certain frequency ranges, known as bandgaps, through a phenomenon known as Bragg scattering. These bandgaps are constrained by the unit cell dimensions, only occurring at low frequencies with wavelengths of the order of the lattice constant [1]. Thus, the application of basic metamaterials was limited to large structures. To expand the use of metamaterials to smaller components, Liu et al. introduced local resonators, showing that locally resonant metamaterials are able to control vibrations at wavelengths much smaller than the lattice constant [13]. With local resonators, bandgaps can result from a combination of Bragg scattering and local resonance. In this case, the resonator parameters have the larger effect on bandgap formation [14]. Indeed, it is also possible for bandgaps to be formed exclusively from local resonances, without the need for Bragg scattering [15]. In addition, vibration control can be achieved at multiple frequencies by adding multiple resonators with different local resonances [16,17].
Bandgaps can also be introduced by the incorporation of piezo-electric materials and shunt circuits. Piezo-electric materials have long been utilized in both active and passive methods for vibration control [18–23] and metamaterials [24–26]. By including piezo-electric elements in a metamaterial, the mechanical system dynamics can be coupled to an easily modifiable shunt circuit. Because the shunt circuit directly affects the structural properties of the coupled piezo-electric material [27], this enables convenient adjustment of the metamaterial's overall properties, enabling optimization of vibration control and energy harvesting, as shown by Abdelmoula and Abdelfiki [28]. Incorporating piezo-electric materials and shunt circuits enables techniques such as the use of negative capacitance [29,30], synchronized switch damping [22,23], or resonant shunt circuits [31,32] to control vibrations and create, tune, or broaden bandgaps [28,33,34]. In addition, shunt circuits also offer an avenue for simultaneous energy harvesting. Tang and Yang investigated the energy harvesting performance of such a piezo-electric energy harvesting model [35], while Hu et al. studied the vibration suppression of an acoustic-elastic metamaterial with embedded piezo-electric element [34].
Combining the two previously discussed methods, researchers have also investigated metamaterials with local resonators coupled to shunt circuits. Some circuits, such as inductor capacitor and resistor-inductor-capacitor (RLC) circuits, can possess resonance of their own, adding to the mechanical resonance of the local resonator itself. Li et al. experimentally confirmed the potential for energy harvesting and simultaneous mechanical wave filtering in a metastructure with local resonators and piezo-electric elements [36]. Hu et al. analyzed a metamaterial beam with piezo-electric resonators using analytical and finite element methods, demonstrating advantages in both vibration attenuation and energy harvesting [37]. Meanwhile, Sugino et al. investigated locally resonant piezo-electric metastructures with a focus on bandgap control [38]. An important parameter in these works and others involving piezo-electric materials is the system's electromechanical (EM) coupling factor. This parameter is dependent on the piezo-electric coupling coefficient, which is itself dependent on the design and material properties of the piezo-electric component [38]. Though this piezo-electric coupling coefficient is usually of the order of for engineering applications [39], signifying weak electromechanical coupling, some features may only be apparent in the case of strong electromechanical coupling. For instance, Hu et al. showed that a system with electromechanical local resonator can manifest a double-valley phenomenon, adding a tunable bandgap to the band structure at high coupling values [34].
This paper applies analytical and numerical methods to investigate the effect of electromechanical coupling and shunt circuit parameters on wave propagation and energy harvesting in a nonlinear acoustic metamaterial with resistance-inductance shunt. Unlike the existing works described above, the proposed metamaterial incorporates nonlinearity in the form of nonlinear springs between the metamaterial cells. This design is similar to a nonlinear metamaterial studied by Bukhari [40]; however, it differs in the design of the shunt circuit. While Bukhari's design included only a shunted resistor, resulting in a resistor-capacitor circuit, the proposed metamaterial includes both resistor and inductor in the shunt circuit, resulting in an RLC circuit and leading to significant changes in the system behavior. This work models the proposed metamaterial through a nonlinear system of governing equations, and the perturbation method of multiple scales is utilized to derive an approximate solution. This solution is validated against the direct analytical solution for the linear case, and the directly integrated numerical results for the nonlinear case. Dispersion relations acquired from the method of multiple scales show that the band structure of the proposed metamaterial significantly differs from the resistance-only metamaterial. The effects of strong coupling and shunt circuit parameters are further explored through parametric analysis of the linear and nonlinear band structures, supplemented by the imaginary band structure for the linear system and numerical frequency response functions (FRFs) for the nonlinear system. Comparisons are made to relevant works that examine RLC shunt systems using shunt theory, showing good agreement with the current results. Further observations not covered by shunt theory are also noted. Following this, the system is excited and numerically simulated to further study the effect of shunt circuit parameters on the voltage response of the metamaterial. Spectro-spatial analysis is applied to the resulting data to explore nonlinear phenomena and energy harvesting potential. Investigation is focused on the effect of shunt parameters, nonlinearity, and the resultant applications for simultaneous vibration control and energy harvesting.
2 Mathematical Modeling of the System
The system studied in this work is a nonlinear acoustic metamaterial coupled with electromechanical local resonators. The metamaterial consists of a chain of cells connected by nonlinear springs as shown in Fig. 1. Each nonlinear spring has linear spring coefficient and nonlinear spring coefficient Λ. Each resonator consists of a substrate covered by a piezo-electric layer, with total effective mass and effective linear stiffness The piezo-electric layer is shunted to an RLC circuit as shown in Fig. 1. This circuit has voltage difference , resistance , and inductance . The piezo-electric layer has capacitance and EM coupling coefficient , as defined by Erturk and Inman in [39]. The absolute displacement of cell n is , and the absolute displacement of the attached piezo-electric resonator is .
A tentative design for the proposed metamaterial chain is shown in Fig. 2. The main cells, modeled in Fig. 1 as point masses with mass , are connected by thin beams to achieve cubic nonlinearity due to midplane stretching as in Refs. [41,42]. A similar design was proposed and experimentally validated by Mojahed et al. in Ref. [43]. In this case, it was shown that a wire of 1080 spring steel with thickness of the order of m (0.006 in) resulted in a nonlinear stiffness of the order of , demonstrating the ability to induce significant nonlinearity through such clamped wires. While the current work is concerned with weak nonlinearity, it has been established that for small deflections, the coupling wire behaves as a Euler–Bernoulli beam, and thus the linear and nonlinear stiffness can be manipulated through the thickness-to-length ratio [43]. Other methods for systems with nonlinear connections that have been experimentally validated include periodic strings [44] and perforated sheets [45].

Tentative design for nonlinear metamaterial cells and resonators: (a) chain and (b) resonator cross section
where and are the piezo-electric constant and permittivity constant, respectively; is the distance from the center of the beam to the center of the piezo-electric layer, is the mass-normalized eigenfunction of the rth vibration mode; and and are the mass per unit length and bending stiffness of the bimorph beam, respectively. Given the material properties of the piezo-electric and substrate layers, these expressions allow calculation of the resonator parameters for any set of bimorph dimensions. These expressions were used to explore the feasibility of the above parameter values.
For preliminary calculations, the piezo-electric layer was taken to be PZT-5A, and the substructure to be aluminum. Equations (1)–(4) were then used to determine the dimensions necessary for desired parameters. For instance, dimensions with length of the order of m, of the order of m, and of the order of m results in weak coupling of the order of N/V and of the order of F. Beam dimensions can be further adjusted to vary selected parameters for future analysis. For example, from Eqs. (1) and (2), it is clear that increasing while decreasing can increase while maintaining the value of . Mass and stiffness values of the resonator will also be affected; thus, the size, material, and design of the main cell may need to be adjusted to maintain the desired mass and stiffness ratio. Additionally, to achieve strong coupling of the order of N/V, resonator dimensions included length of the order of m and of the order of m. It should also be noted that the above capacitance is the modal capacitance, corrected from the clamped capacitance used in Ref. [39] through use of a correction factor to account for the influence of nonresonant modes. This factor is vital for accurate calibration and tuning of the resonator system, and readers are encouraged to consult the more detailed formulation in the work of Høgsberg and Krenk for a full expression [46]. However, as that work shows, is most often of order or less, and thus does not affect the order of magnitude considered by these preliminary calculations.
where is the relative displacement of piezo-electric local resonator n with respect to cell n.
is the electric resonance frequency.
where and
3 Analytical Solution
3.1 Linear Problem.
Solving Eq. (23) reveals six roots for in the form of three complex conjugate pairs. Consequently, the band structure may have up to three pass bands depending on the system parameters.
3.2 Nonlinear Problem.
where .
where denotes the complex conjugate of and .
where and are constants dependent only on system parameters.
3.3 Approximate Solution of Slow Flow Equations.
The nonlinear slow flow equations Eqs. (39) and (40) must be solved to develop the approximate solution for A(T1) and the correction factor for the nonlinear frequency. However, the solution of these equations depends on the values of and . Therefore, we plot their values in Fig. 3 over a frequency range that includes all pass band frequencies. The cases of weak and strong EM coupling, N/V and N/V, respectively, are both considered. Shunt parameter values are Ω, H, and F. The mass of each main cell is kg, and the linear spring constant is N/m. The mass ratio between each piezeoelectric resonator and main cell is , and the mechanical natural frequency of the resonator, , is tuned to the main cell natural frequency such that rad/s. These parameters are chosen based on the similar system studied by Abdelmoula in Ref. [28] and supported by preliminary calculations described above in Sec. 2.
However, when examining the case of strong EM coupling, N/V, it can be seen in Fig. 3(b) that while the magnitude of is less than that of , it is not negligible in all regions; unlike the resistance-only system.
where . The value under the square root, and thus the validity of Eq. (44) as a solution for the system, depends on the value of . If it is negative, then the value under the square root is always positive, so has a valid solution for the system. If is positive, then is only valid for or the trivial solution.
To assess the stability of possible solutions for , the phase portrait of and is plotted for negative and positive values of in Fig. 4. Here, it can be seen that the system goes to zero when is negative and departs from zero when is positive. However, the latter is only defined when . From this, we conclude that the system's only stable solution is the trivial solution of , resulting in Eq. (43) as described above.
The linear and nonlinear dispersion curves derived by the method of multiple scales are then validated against a numerical solution and used to conduct band structure analysis of the system.
4 Band Structure Analysis
4.1 Validation.
To validate the linear and nonlinear dispersion relations, they are compared to a numerical simulation of the system. The system is as described in Fig. 1, with 500 cells coupled to electromechanical resonators with RL shunt circuits. Cells are connected with both linear and nonlinear springs. Parameter values are the same as those used in Sec. 3.3, with weak EM coupling N/V.
where n is the cell number, is the Heaviside function, and is the number of cycles, set as . Excitation amplitude was set to .
For a given wavenumber, the system is simulated for a long time to allow the wave to propagate through the chain. After this, a 2D fast Fourier transform is applied to the time data collected in the wavenumber and frequency domains. The transformed data are then used to determine the natural frequency of the system by finding the frequency associated with the maximum power density point. By determining the natural frequencies corresponding to a range of wavenumbers over the first Brillouin zone, the dispersion curves are numerically constructed and compared to the analytical solution in Fig. 5, which plots normalized frequency ω against wavenumber .
Both linear and nonlinear dispersion relations are compared in Fig. 5. The problem is made linear by setting the nonlinear parameter , resulting in the comparison in Fig. 5(a). For the nonlinear comparison, the nonlinear parameter is set to , resulting in the comparison in Fig. 5(b). These comparisons show good agreement between analytical and numerical solutions for both the linear and nonlinear dispersion relations. It should be noted that the numerical solution is unable to capture certain frequencies of the nonlinear solution in the medium wavenumber region. These areas, known as pseudo-bandgaps, are due to a significant frequency shift associated with transient wavepacket excitation [49]. When excited within these pseudo-bandgaps, the solution appears instead at frequencies within the long and short wavelength limits. This frequency shift will be explored in more detail in Secs. 4 and 5.
4.2 Linear Band Structure.
Next, the validated analytical dispersion relations are used to study the effect of selected parameters on the linear band structure using dispersion curves calculated through Eq. (23). Band structures for both real and imaginary wavenumbers are shown and examined. The parameters examined are EM coupling coefficient and shunt parameters, including resistance , piezo-electric capacitance , and inductance . Unless otherwise noted, parameter values are the same as in Sec. 3.3. In addition, comparisons are also made to the band structure of the system examined by Bukhari et al. in Ref. [40], which is similar except for the absence of inductor in the shunt circuit. Another useful point of reference is the theory of resonant (RL) shunts as described by Thomas et al. [50] and other works in the literature [51,52], which provide support and background for many of the following results.
Varied EM coupling is shown in Fig. 6, with the value of ranging from weak EM coupling at N/V to strong EM coupling at N/V. The previous work by Bukhari et al. Defines “weak coupling” as N/V and “strong coupling” as N/V [40]. The results in Fig. 6 support this definition, as the band structure analysis shows a clear shift as coupling is increased above N/V. For weak coupling values, N/V, a single bandgap forms at . This results in acoustic and optical mode branches similar to the case of resistance-only shunt circuit. However, at N/V, a third mode branch can be observed forming between the acoustic and optical modes, with two smaller bandgaps instead of the single bandgap of the resistance-only system. This split is also visible in the imaginary band structure displayed in Fig. 6(b). As established by the literature on local resonators, these bandgaps correspond to the resonant frequencies of the local resonator [53]. The presence of two resonant frequencies and is due to the coupled RL shunt circuit, which can add a second resonant frequency depending on parameter values such as [50]. Figure 7(a) shows the effect of EM coupling on . They are very close for low coupling values, resulting in a single bandgap. As EM coupling increases, the resonator natural frequencies begin to diverge around N/V, and the bandgaps move with them. The diverging of the resonator frequencies is also visible in Fig. 6(b)'s imaginary band structure, which also gives insight into the effects of on attenuation in the bandgap. Attenuation is largest for weak coupling N/V, decreasing significantly for strong coupling N/V as the bandgap begins to split. However, as is further increased to N/V, the attenuation in the upper bandgap increases while attenuation in the lower bandgap decreases. To further investigate this three-mode system, strong coupling N/V is used in the following analysis unless otherwise noted.

Effect of system parameters on EM resonator coupled natural frequencies, , , H, and F: (a) EM coupling, Ω and (b) shunt resistance, N/V
To determine the effects of shunt resistance on the three-band structure, real and imaginary band structures for varying are shown in Figs. 8(a) and 8(b), respectively. Through these figures, it is clear that shunt resistance affects the formation of the central mode branch as well as attenuation within the bandgaps. At low resistance Ω, two bandgaps are present, while for Ω, this reduces to one. In addition, Fig. 7(b) shows that the lone bandgap is not located exactly at , unlike the resistance-only shunt studied in Ref. [40]. This offset can also be observed in the theory of RL shunts. For low or zero resistance, it has been shown that two bandgaps will be present, located at the coupled EM resonator frequencies . As resistance increases, these two bandgaps will merge into one. The frequency offset of this open-circuit bandgap is dependent on the EM coupling value , indicating that it can be tuned by adjusting this parameter. However, shunt circuit parameters such as and will not affect the location of this bandgap [50]. The imaginary band structure in Fig. 8(b) shows that shunt resistance also affects attenuation within the bandgaps. While high attenuation is beneficial for vibration reduction, this could also drastically reduce energy harvesting due to power harvested being dependent on resonator displacement. Shunt resistance also affects the potential for energy harvesting, due to power harvested being proportional to resistance.

Effect of shunt resistance on linear band structure, , , N/V, H, and F: (a) real and (b) imaginary
Next, the effective piezo-electric capacitance was also varied while keeping all other parameters fixed. Changing capacitance affected the electrical resonance frequency as well as the coupled natural frequencies . The effect of on the latter is illustrated in Fig. 9(a), which shows that only capacitance values F can result in two separate bandgaps. For F, only a single bandgap at is present. Further investigation was carried out by varying relative to the default F, resulting in the real and imaginary linear band structures shown in Figs. 10(a) and 10(b), respectively. This range was chosen to stay close to established values of for piezo-electric elements and produce two distinct bandgaps in the examined nondimensional frequency range of . The real band structure confirms the presence of two bandgaps and three mode branches, and demonstrates how affects the locations of the bandgaps. As increases and decreases, both bandgaps shift to lower frequency. is also shown to influence the width of both bandgaps, with lower values of producing broader bandgaps. The imaginary band structure reveals more about the effects of capacitance on attenuation. As seen in Fig. 10(b), as increases, the attenuation in the upper bandgap increases while the attenuation in the lower bandgap decreases. For the values of and , the band structures show very small lower bandgaps while attenuation increases within the upper bandgap. Meanwhile, smaller values of , such as , show that it is possible to achieve significant attenuation in both bandgaps. While this analysis demonstrates the potential for tuning resonator frequencies and bandgaps using alone, it is important to note that are also affected by , which is fixed at H for these cases. Altering the value of L may affect both the band structure of the system and attenuation within the bandgaps.

Effect of system parameters on EM resonator coupled natural frequencies, , , N/V, and Ω: (a) piezo-electric capacitance, with H and (b) shunt inductance, with F

Effect of piezo-electric capacitance ( F) on linear band structure, , , N/V, Ω, and H: (a) real and (b) imaginary
To investigate the effects of alone, was varied around the default value of H while . These values of were selected based on Fig. 9(b), which shows that shunt inductance H is required for two distinct bandgaps in the studied frequency range. For smaller , only a single bandgap at is present. The resulting real and imaginary band structures are shown in Fig. 11 and match the trends shown in Fig. 9(b), with bandgaps growing farther apart as inductance increases. Similar to the effects of , broader bandgaps are present for smaller values of . In addition, both bandgaps shift to lower frequency as increases, also similar to the effect of . However, the locations of the band gaps are not the same for the same change in and . Another difference between the parameters' effects can be seen in Fig. 11(b). As is decreased, attenuation in the upper bandgap decreases. However, as is decreased, attenuation in the upper bandgap increases. This shows that while the bandgap locations may be adjusted by either or , improving performance of the system may require further optimization of one or both parameters.

Effect of shunt inductance ( H) on linear band structure, , , N/V, Ω, and = 1.13 × 10−10 F: (a) real and (b) imaginary
Finally, both and were varied to investigate their combined effect on the system. Inductance and capacitance are chosen to keep the electrical resonance frequency fixed at . Since , any change in requires an opposite change in to keep constant. The resulting band structure is displayed in Fig. 12. The effect of keeping fixed is readily apparent, with several significant differences from the previous analyses. Figure 12(a) shows that all band structures have three propagation modes, with the second mode centered around ω = 1. As L increases, the upper and lower bandgaps shift to higher and lower frequency, respectively, with the upper bandgap moving farther than the lower. As seen in Fig. 12(b), the upper bandgaps have higher attenuation than the lower bandgaps, but attenuation for both bandgaps increases with increasing L. This is in contrast to the previous analyses, where it was found that shifting L or resulted in one bandgap gaining attenuation while the other lost attenuation. However, the width of the bandgaps is similar to this pattern, with increasing L causing the upper bandgap to broaden while the lower bandgap grows more narrow.

Effect of shunt inductance ( H) on linear band structure with tuned resonator , , N/V, and Ω: (a) real and (b) imaginary
4.3 Nonlinear Band Structure.
Following the investigation of the impact of shunt circuit parameters on the linear band structure, the nonlinear band structure was also studied by taking into account the correction factor defined by Eq. (42). The effects of the nonlinear stiffness parameter are varied to model chains with nonlinear hardening () and nonlinear softening (). The effects of nonlinear hardening and softening are then combined with the effects of EM coupling and other shunt parameters.
The effect of nonlinear hardening or softening on the dispersion curves with default parameters is shown in Fig. 13(a). As expected, hardening shifts portions of the band structure to higher frequency, while softening shifts them to lower frequency. The magnitude of the shift is negligible for small wavenumber (long wavelength limit) as well as low frequencies. It is more pronounced for the third mode branch, especially at larger wavenumbers (short wavelength limit). This is similar to the resistance-only system, which showed larger effects of nonlinearity within the optical mode [40]. However, there is a difference when examining the lower two modes. For the resistance-only system, effects of nonlinearity were greatest at large wavenumbers for both branches. However, in the RL metamaterial, the effects of nonlinearity on the first and second modes decrease near the bandgaps, resulting in the effects of nonlinearity being largest in the medium wavelength region. The cases of linear and nonlinear chains with weak EM coupling are displayed in Fig. 13(b). Here, the effects of nonlinear hardening and softening on the mode cutoff frequencies are largely the same as the case of strong coupling. The band structure has two branches rather than three, matching the effects of weak coupling previously established, and being similar to the resistance-only system. Like the strong coupling case, the first mode differs from the resistance-only system by having low effects of nonlinearity at high wavenumber. Overall, the effects of nonlinearity are most prominent in the third mode at high wavenumber and frequency.

Effect of system parameters on nonlinear band structure, , , and: (a) default parameters, N/V, , , (b) weak EM coupling N/V, (c) , untuned resonator , and (d) , tuned resonator
Next, the shunt inductance is increased to in Fig. 13(c) while keeping . This case corresponds to one of the band structures in Fig. 11. There, two bandgaps were shown, with the upper bandgap approaching and the lower bandgap around . These are also visible here. Similar to the linear case, the lower bandgap near is very narrow. It has no significant interaction with the effects of nonlinearity, which is expected due to its low frequency/wavenumber value as well as its small bandwidth. The upper bandgap is similar to the single bandgap of the low coupling case in Fig. 13(b), though it is slightly broader, and shifted to higher frequency. The effects of hardening and softening nonlinearity are also similar to the low coupling case for the upper modes. The case of increased shunt inductance and tuned resonator, , is shown in Fig. 13(d). As seen in Fig. 12, increasing while keeping the resonator tuned results in two bandgaps and three modes, with the second mode branch covering an increased frequency range compared to the default parameter case in Fig. 13(a). Due to this, the effects of nonlinearity on the second mode are more prominent for the increased inductance case, especially in the medium wavelength limit. However, effects of nonlinearity still decrease near the bandgaps, so nonlinearity does not significantly affect the location of the bandgap boundaries. Like the previous cases, the effects of nonlinearity are most pronounced in the large wavenumber region of the third mode.
where s. The resulting FRFs, shown in Fig. 14, support the observations drawn from the band structure analysis while also making other details more apparent. Notably, the frequency response significantly differs from that of a corresponding single-cell structure, such as those studied in shunt theory. While a single-cell structure would have peaks, with high transmissibility localized to a number of individual frequencies, the periodic system studied here displays pass bands, with high transmissibility present over a number of frequency ranges. This is of great benefit for both vibration control and energy harvesting. This observation is also supported by the imaginary band structures, which all show flat attenuation values between the bandgaps.
Two cases are included here, both with tuned resonator, . Figure 14(a) shows the default case, with and . Figure 14(b) shows a case of increased inductance, with and . For both cases, the numerical FRF concurs with the analytical band structure in Fig. 13, showing bandgaps at the same frequency locations. The effects of nonlinearity are also visible on the upper bound of the third mode, with a shift in cutoff frequency caused by nonlinearity. The chain with nonlinear hardening shifts the upper bound of the mode to higher frequency while softening shifts to lower frequency. This confirms the effects shown in the analytical band structure. Furthermore, the FRF provides insight into wave propagation within the second mode not visible through band structure examination or shunt theory. When , Fig. 14(a) shows visible attenuation in the second mode, with transmissibility several orders of magnitude below that in the other two modes. In Fig. 14(b), some minor attenuation is also visible, but the transmissibility in the second mode is closer to that in the other modes. Thus, it can be concluded that decreasing L will result in increased attenuation within the second mode itself. This behavior is not revealed by any of the previous analysis, and adds another point of interest to the three-mode band structure.
5 Wavepacket Excitation Analysis
To further study the system response within the previously described band structure, the linear and nonlinear chains are excited by the transient wavepacket defined in Eq. (45). As this wavepacket is a function of wavenumber, the reciprocal of wavelength, this enables us to examine the system in multiple wavelength regions. The selected wavelength regions are the long wavelength region (), the medium wavelength region (), and the short wavelength region (). Each wavelength region contains multiple mode branches with frequency–wavenumber relations shown in the preceding band structure analysis. Thus, we can examine each mode in a selected wavelength region, determining the behavior of the system for varying wavenumber and the corresponding frequency values. It should be noted that the linear dispersion relations were used as an approximation to relate the wavenumber to the corresponding frequency in each propagation mode [49].
5.1 Spatial Profile.
First, the system is simulated for 1500 s (nondimensional time) to allow the wavepacket to propagate through the chain. After this period, the spatial profile of the chain's voltage response is studied to examine the effects of nonlinearity. The system is excited within all three modes, and multiple wavenumber regions. Following from the previous analysis, two cases are studied: and , both with tuned resonator. The frequency–wavenumber relations for these cases can be found in the corresponding band structures in Figs. 13(a) and 13(d), respectively. Selected voltage response profiles for linear and nonlinear chains are shown in Figs. 15–17, with the input (initial excitation) also included for reference. All voltage values are normalized to the peak input voltage.

Spatial profile of voltage response, first mode branch, for , N/V, , Ω: (a) , , (b) , , (c) , , and (d) ,

Spatial profile of voltage response, third mode branch, for , N/V, , Ω: (a) , , (b) , , (c) , , and (d) ,
Figure 15 shows the response of chains excited within the first mode. Here, the system is excited in the long wavelength region (); as well as the medium wavelength region (k = π/3), where the mode approaches the bandgap. In the long wavelength region, Fig. 15(a) and 15(b), both cases show similar results, with effects of hardening and softening nonlinearity visible but not significant. This small effect of nonlinearity is consistent with the band structure shown in Fig. 13. All values of inductance and nonlinearity show vibration attenuation, with the output voltage profiles an order of magnitude less than the input. In the medium wavelength region, the voltage response is significantly reduced for both cases, indicating that vibration attenuation increases as the bandgap is approached. This is supported by the similar observations from the linear band structure examined in Fig. 12. Again, effects of nonlinearity are visible, but not significant. Comparing Figs. 15(c) and 15(d), the case of increased inductance shows higher voltage amplitude, approximately double of the default case. However, both are two orders of magnitude smaller than the input. When compared to the corresponding response of the metamaterial with resistance-only shunt and low EM coupling studied by Bukhari, there is a clear difference in the effects of nonlinearity [40]. In that system, there is negligible effect of nonlinearity in the long wavelength region, but significant effects in the medium and short wavelength regions, as opposed to the low effects in all regions for the current system. In addition, the low-coupling system does not show the decrease of vibration propagation in the short wavelength region adjacent to the bandgap.
The response of the second mode is shown in Fig. 16. For this mode, only the medium wavelength region () is shown here. However, the long and short wavelength regions were examined and found to show strong vibration attenuation and low effects of nonlinearity consistent with the observations made for the first mode in Figs. 15(c) and 15(d). At , both the default and increased inductance cases have output two orders of magnitude less than the input. Since the medium wavelength region is not near either edge of the pass band, the low wave propagation is not an example of increased attenuation near the bandgaps. It can be concluded that even though the band structure in Fig. 13 shows a pass band, the second mode branch has low vibration propagation throughout its frequency range. This is consistent with the reduced transmissibility for the second mode observed in Fig. 14(b). As in the medium wavelength region of the first mode, the increased inductance case has roughly doubled output magnitude, corresponding to the increased transmissibility in Fig. 14(b). Despite the low overall amplitude, this mode shows significant effects of nonlinearity for both cases. The wave profile of the chain with nonlinear hardening shows the presence of narrow solitary waves, more localized than the linear response and with slightly greater amplitude. The softening chain response is more dispersive than the linear and hardening chains, with a lower amplitude spread out across more of the chain. Overall, in the medium wavelength region, effects of nonlinearity are similar to those observed for the second mode in the resistance-only system [40]. However, the response of the resistance-only system was much larger, matching or exceeding the input amplitude. In addition, the resistance-only system showed larger effects of nonlinearity in the short wavelength region (), in contrast to the small effects in the current case.
In Fig. 17, spatial profiles of the response to excitation within the third mode are shown. The response in the long wavelength region () is not included, but shows low amplitude and effects of nonlinearity near the lower bound of the pass band, similar to Figs. 15(c) and 15(d). Figures 17(a) and 17(b) show the response of the two cases to excitation in the medium wavelength region. Output amplitude is larger than the second mode, though still not as large as the input. The increased inductance case shows slightly smaller output. Effects of nonlinearity are similar to the medium wavelength region of the second mode, with hardening resulting in higher-amplitude solitary waves and softening resulting in lower-amplitude dispersive components. In the short wavelength region (), the output amplitude increases, with even the linear and softening responses having amplitude similar to the input. In both cases, the hardening chain response contains a solitary wave with amplitude larger than the input. Again, the increased inductance case has smaller output than the default parameter case. The wave profiles in the third mode are most similar to the response of the second mode in the resistance-only system [40]. While the medium wavelength region still shows smaller response amplitude compared to the resistance-only system, the response in the short wavelength region includes solitary waves with even larger amplitude than the corresponding ones in the resistance-only system. It should be noted that the previously described vibration attenuation near the bandgaps is not present at the upper bound of the third mode. Thus, it becomes clear that this attenuation is only present near the bandgaps between the modes; i.e., the bandgaps resulting from the presence of the shunt inductor. Combined with the previously described vibration attenuation across the second mode, this results in a large range with either bandgaps or very low wave propagation, an observation with great potential for vibration control.
5.2 Spectro-Spatial Analysis.
Another phenomenon that merits study is the frequency shift noted in Sec. 4.1. In Fig. 5, it was observed that when excited by a transient wavepacket and solved numerically, the medium wavenumber region of the upper mode branch does not appear in the predicted band structure. This indicates a significant frequency shift in that region, supported by similar observations by Zhou [49] and Bukhari [54]. To investigate possible effects of shunt parameters on these frequency shifts and the energy harvesting potential, spectro-spatial analysis is utilized by applying signal processing techniques to the results of the previous wavepacket excitation. We apply a short-term Fourier transform (STFT) to the data in the spatial-temporal domain in order to study any effects that manifest in the spatial wavenumber domain. The resulting data are presented in spectrograms with cell number on the x-axis and wavenumber on the y-axis. The color of the spectrograms represents the relative voltage amplitude, and both input and output (harvested) signals are displayed on each figure for purposes of comparison. A Hann window was used for this STFT, with window length matching the size of the input burst. Window length was defined as equal to the wavelength () of the input wavepacket. This window size was chosen in order to contain the short spatial components and for convenient monitoring of the wave distortion [49,55]. This window is represented by the dashed lines in each figure. Similarly to the spatial profile analysis, the system is examined for selected modes within a wavelength region. Here, we display the results for the second and third modes in the medium wavelength region in order to investigate properties such as frequency shift. The frequency–wavenumber relations for these modes can be seen in the band structures shown in Figs. 13(a) and 13(d) for and , respectively.
The results for the second mode branch in the medium wavelength region are shown for both cases in Fig. 18. For the default parameter case, the voltage input and output of the linear chain are displayed in Fig. 18(a). It is clear that the voltage output differs slightly from corresponding results of the system studied by Bukhari [54]. This is due to the addition of inductance to the shunt circuit as well as the increased EM coupling. The output voltage profile shows significantly less energy than the input, an observation consistent with the previous analysis of the spatial profiles. However, despite the low energy, the results for hardening and softening chains in Figs. 18(c) and 18(e) show frequency shifts similar to those observed by Bukhari [54]. Specifically, the outputs for hardening and softening chains show different frequency distributions of their energy. For the hardening case, the high-energy solitary component near cell 200 shifts to higher frequency, close to the border of the Hann window. Two dispersive components with lesser energy are also present at the middle of the window, along with an additional low-energy component on the edge of the window. For the softening chain, two dispersive components are present, one inside the Hann window and one outside. Overall, the response is more dispersive than Bukhari's system without inductor and has lower output voltage [54]. When shunt inductance is increased to , the output voltage profiles become more dispersive for all three chains, and gain more energy. Figure 19(b) shows energy harvested for the linear chain, with a voltage response over half the chain rather than the more localized response of the default parameter case. This dispersive component shows an increase in the energy harvested, though still not to the level of the input. For the hardening chain, shown in Fig. 19(d), the response includes one frequency-shifted solitary wave as in the default parameter case. This localized solitary wave, located near cell 300, has a larger frequency shift than the default case in Fig. 19(c), with half the energy outside the Hann window. In addition, two other components between cells 400 and 500 have become more localized than the default case and have higher energy, similar to the shifted solitary wave. For the softening chain in Fig. 19(f), the energy components are overall more dispersive than the default case, with most of the energy within the Hann window. These results are a good match to the spatial profile in Fig. 16, confirming the higher energy harvested by increasing inductance and showing that the frequency shift observed in similar systems is still present despite the overall decrease in wave propagation.

Short term Fourier transform of voltage response, second mode branch, , for , N/V, , Ω: (a) linear chain, , (b) linear chain, , (c) hardening chain, , (d) hardening chain, , (e) softening chain, , and (f) softening chain,

Short term Fourier transform of voltage response, third mode branch, , for , N/V, , Ω: (a) linear chain, , (b) linear chain, , (c) hardening chain, , (d) hardening chain, , (e) Softening chain, , and (f) softening chain,
Figure 19 displays the voltage profiles when both cases are excited in the medium wavelength region of the third mode branch. The majority of the observations made for the second mode still hold true here, with the effects of hardening and softening nonlinearity still present, including frequency shift. However, in this wavenumber/frequency region, more energy is harvested in the default case, rather than the increased inductance case. This corresponds to the spatial profile results in Fig. 17.
6 Discussion
Here, we summarize and discuss the results of this work. First, linear band structure analysis explored the effects of the RL proposed shunt circuit. Parametric analysis was carried out by varying EM coupling , shunt resistance , piezo-electric capacitance and inductance ; and observing the effect on the real and imaginary band structures for the linear system. This study showed that there are three main options for bandgap placement: a single bandgap at for weak EM coupling; a single bandgap offset from for strong EM coupling and high resistance; and two bandgaps for strong EM coupling and low resistance. For each of these cases, the location and attenuation of these bandgaps can be tuned using certain resonator parameters. In addition, it was shown that weak EM coupling causes the system's behavior to be similar to a metamaterial with resistance-only shunt. Overall, these observations of the linear system were supported by the literature on shunt theory.
Using the derived nonlinear solution and numerical methods, analysis then moved forward, providing insight not available through linear analysis or shunt theory. While the cubic spring nonlinearity within the proposed metamaterial has little direct interaction with the EM resonator and shunt, the system's behavior includes effects of both nonlinearity and the resonator, as seen in the nonlinear band structure analysis. In the dispersion curves of the three-mode system, the effects of nonlinearity were most prominent in the third mode branch at the short wavelength limit, similar to systems with resistor-only shunt previous studied in the literature. However, for the lower two modes, the effects of nonlinearity are most significant in the medium wavelength region, and decrease near the inner bandgap boundaries. In addition, study of transmissibility FRF diagrams for the nonlinear system also revealed that depending on parameter values, vibrations may be attenuated around the inner bandgap boundaries and within the central mode. These FRFs also emphasized the advantage of a periodic system unlike those studied in simple shunt theory.
Through numerical simulation of the system, the responses of linear and nonlinear chains were further explored using wavepacket excitation and spectro-spatial analysis. The spatial profiles of the voltage response supported previous observations, including the low vibration propagation and low influence of nonlinearity in certain areas. Specifically, these effects were found to be focused around the inner bandgaps. Wave propagation was reduced in the portions of all three mode branches adjacent to these inner bandgaps. The second mode, located between these bandgaps, was shown to have strong vibration attenuation across the pass band. In addition, the effects of nonlinearity decreased for all three mode branches as they approached these inner bandgaps at both small and large wavenumber; marking a difference from the similar resistance-only system, which showed the effects of nonlinearity always increasing at large wavenumber. The response for the lower two modes was found to be less than the corresponding results for the resistance-only system. Increasing shunt inductance was shown to result in higher voltage response, though still significantly less than the resistance shunt system. However, the energy harvested in the third mode was higher than the corresponding results for the resistance-only system. In addition, the frequency shift due to nonlinear chain was still observed in the medium wavelength region for multiple inductance values.
In addition to the observations described above, this work also poses questions for future study. The parametric analysis described above provides a basic understanding of the individual effects of shunt circuit parameters, but also shows the presence of more complex relations. For instance, and both influence the band structure, but the effects of each are also dependent on the other's value. While this work adjusts these parameters both apart and together, further study is required to optimize their values following shunt theory. In addition, while this paper includes a preliminary design of the proposed chain and general calculations for achieving desired parameter values, this design will be refined in future work. This will lead to the physical modeling of the proposed nonlinear metamaterial for experimental work.
7 Conclusions
This work contains the analysis of a nonlinear, electromechanical metamaterial coupled to a shunt circuit with both resistor and inductor. The system consists of a chain of cells connected by nonlinear springs, with each cell coupled to an electromechanical resonator consisting of piezo-electric element and shunt circuit. Both the cells and resonators were modeled as spring-mass systems, with the resonator system also coupled to the dynamics of the shunt circuit. The resulting governing equations were solved analytically using the perturbation method of multiple scales, and the approximate solution was validated by direct numerical integration. The validated analytical solutions were then used to examine the system's band structure to view the effect of electromechanical coupling and resonator parameters including mass ratio, shunt resistance, and shunt inductance. Observations of the linear chain were found to be consistent with the literature on local resonators and shunt theory. For the nonlinear case, band structure analysis was also supported by transmissibility diagrams obtained by direct numerical solution of the system. It was shown that with certain parameter values, a band structure with three propagation modes can be obtained, rather than the two-mode structure observed for a similar system with resistance-only shunt. This three-mode system was then investigated through parametric analysis, focusing on the effect of shunt parameters on bandgap location as well as vibration attenuation within the pass bands and bandgaps. While the first and third modes of the proposed metamaterial were found to behave similarly to the acoustic and optical modes of the two-branch system, the second mode branch showed potential for vibration attenuation.
The effects of nonlinear springs on the system dynamics, especially the interaction with shunt inductance, were also studied. In addition to the dispersion relations obtained from the method of multiple scales, nonlinearity was also examined through numerical simulation of the metamaterial chain under transient wavepacket excitation, and spectro-spatial analysis of the resulting data. This allowed for analysis of the spatial profile of the linear and nonlinear chain responses, as well as spectrograms of the voltage response. Investigation of the spatial profile confirmed the presence of a region with strong vibration attenuation around the second mode branch and adjacent bandgaps. Voltage response across the first and second modes of the proposed metamaterial was shown to have lower amplitude than the resistance-only system, while the third mode was shown to harvest more voltage than the resistance only system. Increasing shunt inductance increased the energy harvesting potential of both linear and nonlinear chains for the lower two modes, but decreased the voltage harvested in the third mode. Both values of inductance did retain the effects of nonlinearity, including frequency shift. Overall, the addition of an inductor to the shunt circuit resulted in a system with a highly tunable three-branch band structure in addition to the benefits of the previously studied resistance-only system. The added central mode and adjacent bandgaps possess strong potential for vibration control, while the third mode is well suited for energy harvesting.
Funding Data
National Science Foundation (NSF) (Grant Nos. CAREER ECCS1944032 and CMMI-2038187; Funder ID: 10.13039/100000001).
Data Availability Statement
The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.