## Abstract

Energy dissipation in polymeric composite metamaterials requires special mathematical models owing to the viscoelastic nature of their constituents, namely, the polymeric matrix, bonding agent, and local resonators. Unlike traditional composites, viscoelastic metamaterials possess a unique ability to exhibit strong wave attenuation while retaining high stiffness as a result of the “metadamping” phenomenon attributed to local resonances. The objective of this work is to investigate viscoelastic metadamping in one-dimensional multibandgap metamaterials by combining the linear hereditary theory of viscoelasticity with the Floquet-Bloch theory of wave propagation in infinite elastic media. Important distinctions between metamaterial and phononic unit cell models are explained based on the free wave approach with wavenumber-eliminated damping-frequency band structures. The developed model enables viscoelastic metadamping to be investigated by varying two independent relaxation parameters describing the viscoelasticity level in the host structure and the integrated resonators. The dispersion mechanics within high damping regimes and the effects of boundary conditions on the damped response are detailed. The results reveal that in a multiresonator cell, strategic damping placement in the individual resonators plays a profound role in shaping intermediate dispersion branches and dictating the primary and secondary frequency regions of interest, within which attenuation is most required.

## 1 Introduction

The study of periodic systems has recently matured and evolved into a research field that offers a versatile platform for the synthesis of mechanical structures with properties that go beyond traditional and natural limits [1]. Specifically, in the domain of acoustic and elastic wave propagation, periodic architectures such as phononic crystals (PCs) and locally resonant metamaterials (LRMs) have been shown to harbor unique and unprecedented wave dispersion characteristics, which culminate in frequency bandgaps, i.e., extended ranges of forbidden wave propagation, as well as interesting directional filtering and energy guidance capabilities [2–6]. Most recently, the tailoring of such PCs and LRMs, and variations thereof, via active inserts [7,8], topology optimization [9], symmetry breakage [10], fluid-coupling [11], or selective control of material properties in time and space [12] have given rise to unconventional patterns related to insulation [13], atomic and nanoscale heat transfer [14,15], and wave reciprocity [16,17], which uniquely position such periodic systems to fundamentally alter and transform a wide range of vibroacoustic engineering applications.

Although the hallmark features of PCs and LRMs (e.g., frequency bandgaps) take place in the nondissipative form of such systems, the effect of damping cannot be ignored [18,19]. Practically speaking, almost all composite and composite-like structures exhibit a certain degree of viscous or viscoelastic damping, which interferes with and alters their dispersion mechanics and dynamic response. More importantly, owing to natural trends and coupling of mechanical properties (often depicted in Ashby charts of commonly available materials), high dissipation levels are usually associated with softer materials, posing serious questions about their ability to function in certain environments and heavy applications. Over the past decade, LRMs—where the interplay between a host structure and periodic resonant inclusions instigates subwavelength and size-independent bandgaps—have jumped to the forefront of material candidates that can be tailored to exhibit strong attenuation while retaining high stiffness. The ability of LRMs to provide this elusive combination stems from the concept of “damping emergence,” i.e., the ability of an internal structure of a material to whip up its dissipation performance compared to other materials with the same damping composition and equivalent static properties. The latter notion is often referred to as *metadamping*, a term coined by Hussein and Frazier who first conceptualized this problem in canonical spring-mass lattices [20]. The phenomenon has since spurted a number of efforts, which have examined metadamping in the context of viscoelastic metamaterials [21], nontraditional crystals [22], flexural systems [23], as well as fractional calculus [24]. Experimental evidence of positive and negative metadamping (i.e., enhanced and reduced dissipation) in continuous flexural systems was recently presented [25]. A comprehensive overview of metadamping in elastic metamaterials can be found in the study by Bacquet et al. [26].

Fiber-reinforced polymers (FRPs) are a class of composite materials characterized by high stiffness and strength, low density, as well as tunable damping properties that differ from traditional engineering materials. The driving mechanism of energy dissipation in FRPs stems from the polymeric matrix and bonding agents that exhibit viscoelastic features. By increasing the volume fraction of the matrix and/or bonding agents, damping will increase at the expense of stiffness and strength. Thermosets (such as epoxy) are generally preferred over thermoplastics because they exhibit higher stiffness and better adhesion properties balanced with eminent (yet lower) damping characteristics [27,28]. Although considered the main contributor to damping, the matrix is not the only source of energy dissipation in composite materials. Several other sources include [29,30] (1) damping in the fibers, (2) damping due to interphase, (3) damping due to damage, (4) thermoelastic damping, and (5) viscoplastic damping. The current literature contains numerous methods for predicting composite damping, ranging from analytical and computational models [31–34] to experimental characterization techniques [35–38]. When FRP is augmented with viscoelastic resonant inclusions, the resultant viscoelastic metamaterial attains superior attenuation levels combined with a high strength-to-mass ratio. The effects of dissipation in viscoelastic metamaterials have been studied, for example, using Kelvin–Voigt and generalized Maxwell models [39] and in the context of the dynamic homogenization theory [40].

Motivated by the previous expectations and the fact that LRMs can operate efficiently only within stringent damping limits or constraints, this article is concerned with the modeling and analysis of the *metadamping* phenomenon associated with multiresonator viscoelastic metamaterials. To this end, the free wave approach is used to investigate the dispersion mechanics of a graphite-epoxy metamaterial beam both in the infinite and finite medium configurations. This article is organized as follows: following Sec. 1, the finite element formulation for computing viscoelastic band structures is detailed along with numerical validations against commonly used viscously damped models. Section 3, which constitutes the greater part of this effort, focuses on elucidating the metadamping phenomenon encountered in the damping-frequency band structures across independent configurations of one-dimensional (1D) composites comprising two viscoelastic resonators. Within the analyses, the mechanics of dispersion in high damping regimes and the effects of physical boundary conditions on the beam’s damped response are discussed. The obtained results will then be generalized for composite metamaterials with multiple resonators. Finally, the conclusions of the current work are summarized.

## 2 Viscoelastic Band Structure Computation

We start by presenting a comprehensive mathematical framework for computing damped band structures (also known as dispersion diagrams) of viscoelastic and viscous 1D periodic structures based on the free wave approach. Within the analysis, important distinctions between (1) metamaterial and (2) phononic unit cell models are explained in details.

### 2.1 Viscous Damping Model.

**u**. According to Bloch’s theorem, such wave field vector can be expressed as follows:

*a*represents the spatial periodicity of the unit cell in the

*x*-direction,

*κ*is the wavenumber, and

*λ*is a complex frequency function that permits wave attenuation in time. Applying the periodicity constraint, Eq. (1) can be rewritten as a relationship between the unit cell boundaries such that

**u**is discretized into generalized displacements satisfying the motion equations given by $MU\xa8+DU\u02d9+KU=0$, where

**U**is the free generalized displacements vector and

**M**,

**D**, and

**K**denote the unit cell mass, viscous damping, and stiffness matrices, respectively. For discretized models, Eq. (2) is equivalent to a set of equations relating the generalized displacements of the unit cells’ right boundary to their left counterparts. This implies that $U=[ULUIUR]T$, which contains the full set of rearranged degrees-of-freedom (DOFs), can be expressed as $U=PU~$, where $U~=[ULUI]T$ is the periodic generalized displacements vector that contains the minimal set of DOFs (including only internal and left boundary displacements and rotations), and

**P**is the Bloch periodicity matrix, which is expressed as follows:

**I**and

**0**are identity and null matrices of proper sizes, respectively. Consequently, unit cell equations can be written in terms of $U~$ as follows:

*a*)

*b*)

*c*)

*λ*(

*κ*) = −

*ζ*

_{d}(

*κ*)

*ω*

_{r}(

*κ*) ±

*jω*

_{d}(

*κ*), where

*ω*

_{r},

*ω*

_{d}, and

*ζ*

_{d}are the wavenumber-dependent undamped frequency, damped frequency, and damping ratio, respectively. While this research effort is primarily concerned with the analysis of viscoelastic dissipation, the above formulation for viscous wave propagation is included as a benchmark against which the viscoelastic behavior will be evaluated as outlined in Sec. 2.3.

### 2.2 Viscoelastic Damping Model.

**G**(

*t*−

*τ*) is a matrix of damping kernel functions (also known as hereditary functions) with a series of parameters that are typically identified on master curves. Many viscoelastic models have been developed and applied to replicate the frequency-dependent viscoelastic material properties including, but not limited to, the Golla–Hughes–McTavish [42], the anelastic displacement field [43], and the exponential damping models [44]. Some of the kernel functions commonly known in the literature are presented in Table 1. In this article, the exponential damping model is implemented due to its mathematical convenience for modeling complex viscoelastic structures incorporating combination of damping mechanisms arising from the host structure and/or local resonators. The model also enables us to tune the relaxation parameters as needed to replicate the behavior of (a) a purely elastic material, (b) a purely viscous material, or (c) a viscoelastic material that provides a critical validation mechanism of the results presented here.

*μ*

_{k}are known as the relaxation parameters and

*n*is the number of parameters required to replicate the damping behavior of a given viscoelastic material. Substituting Eq. (11) into (10) yields

*a*)

*b*)

*a*)

*b*)

**Φ**

_{kp}and

**Φ**

_{ks}are the eigenvectors corresponding to nonzero (

*p*) and zero (

*s*) eigenvalues of $D~k$, respectively. A reduced set of internal variables $Z~kr$ is then introduced, such that $Z~k=\Phi kpZ~kr$. Upon substituting the previous relation into Eq. (16) and premultiplying Eq. (16

*b*) by $\Phi kpT$, we arrive at

*a*)

*b*)

*a*) and (18

*b*) to be combined into the reduced state-space form given by

*a*)

*b*)

From Eqs. (17) and (20), it can be observed that in the limit $\mu k\u2192\u221e\u2200k$, the viscous state-space form given in Eq. (7) is recovered. The previous implies that high values of *μ*_{k} yield a viscous-like behavior that is less dependent on the history, while low values of *μ*_{k} are more representative of a full-fledged viscoelastic behavior.

### 2.3 Numerical Validations.

Before pursuing an investigation of viscoelastic metadamping, which is a phenomenon associated with dissipative metamaterials, an initial set of results produced from the analysis of viscoelastic phononic materials is briefly outlined here as a baseline for future comparisons. Moreover, this step provides a simple way for validating the viscoelastic model presented in Sec. 2.2. To this end, let us consider a 1D fiber-reinforced composite structure as shown in Fig. 1(a). It consists of a 3 mm thick, 40 mm wide, and 360 mm long [0 deg] graphite-epoxy cantilever beam with the following material properties [49]: *ρ* = 1389.23 kg/m^{3}, *E*_{1} = 144.80 GPa, *E*_{2} = 9.65 GPa, *G*_{12} = *G*_{13} = 4.14 GPa, *G*_{23} = 3.45 GPa, and *ν*_{12} = 0.3, where *ρ*, *E*, *G*, and *ν* denote the density, elastic modulus, shear modulus, and Poisson’s ratio along the different directions, respectively. A series of 12 lumped masses are periodically attached to the lower surface of the composite beam at equally spaced distances *a* = 30 mm, as shown in Fig. 1(b). The lumped masses are chosen such that *m*_{v} = 0.4 *m*_{b}, where *m*_{v} and *m*_{b} denote the lumped and beam masses in a single-unit cell. A notable advantage of the viscoelastic model presented in Sec. 2.2 is its ability to be validated against classical viscous damping models. The previous hypothesis can be validated by assuming that **D**_{b} = *α***K**_{b}, where **D**_{b} and **K**_{b} represent the viscous damping and stiffness matrices, respectively, associated with the composite beam and *α* is a proportionality constant. Throughout this study, the FEM used to obtain unit cell predictions is based on the first shear deformation theory and Timoshenko beam assumptions (see Ref. [50] for details). Furthermore, while fiber orientation plays a significant role in dictating several mechanical properties of a given composite material, its effect on energy dissipation is beyond the scope of this effort. Alternatively, a single relaxation parameter *μ*_{b} will be used to describe the damping behavior of the considered unidirectional composite beam, enabling the metadamping phenomenon to be evaluated in a general sense rather than the behavior of a particular viscoelastic material. As such, it becomes plausible to tune the beam’s viscoelasticity level (dependence on the history) to one of three damping regimes, namely, (1) undamped, (2) viscoelastic, or (3) viscous [21].

In Figs. 2(a) and 2(b), frequency and damping ratio band structures are shown for viscous (*α* = 3.58 × 10^{−6}) and viscoelastic (*μ*_{b} = 2 × 10^{10}) plain composite beams (in the absence of lumped masses). For such a large value of *μ*_{b}, it can be observed that the band structures of both beams are almost identical across the entire wavenumber spectrum. To clearly demonstrate the effects of viscoelasticity on the mechanics of dispersion, Fig. 2(c) is constructed from a combination of Figs. 2(a) and 2(b). The main advantage of the resulting figure (henceforth denoted as the damping-frequency band structure) stems from its readability convenience since it eliminates the wavenumber completely, providing a direct correlation between Bloch damping ratios and the associated oscillatory frequencies. The damping-frequency band structure is computed for three different *μ*_{b} values. The first value is selected to be *μ*_{b} = 2, which is small enough to replicate the behavior of an undamped beam, producing substantially small values of Bloch damping ratios across the entire frequency spectrum as inferred from the flat line coinciding with the frequency axis. For *μ*_{b} = 2 × 10^{10}, which was used to construct Figs. 2(a) and 2(b), the viscous damping case is recovered, hence validating the proposed viscoelastic model for the case of a plain composite beam. The figure also depicts an additional intermediate curve representing a more realistic behavior of a viscoelastic beam with *μ*_{b} = 2 × 10^{5}. For this particular case, modal damping ratios of the corresponding finite cantilever beam are plotted as discrete points in the same figure and showing excellent agreement with the Bloch damping ratios in the considered frequency range (0–20 kHz).

In a second validation step, the same parameters (*α* = 3.58 × 10^{−6} and *μ*_{b} = 2 × 10^{10}) are used to construct the frequency and the damping ratio band structures shown in Figs. 3(a) and 3(b) for the corresponding viscous and viscoelastic phononic composite beams, respectively. Unlike the plain beams, the band structures expectedly capture Bragg bandgaps spanning two different frequency regions (10.14–13.58 kHz and 30.01–36.19 kHz). Again, the viscous and viscoelastic beams produce identical results for a large *μ*_{b} value, and similar observations can be made about Figs. 2(c) and 3(c), with the exception of Bloch damping ratio solutions that are notably absent within the bandgap region in Fig. 3(c). Another important observation is that the existence of a Bragg bandgap in a given phononic beam does not fundamentally influence its damping behavior. Finally, it is worthy noting that “band-overtaking” takes place in both Figs. 2(a) and 3(a), which is a phenomenon associated with higher frequency bands crossing lower frequency ones as a result of increased damping levels [26].

## 3 Metadamping Analysis

The objective of this section is to leverage the aforementioned framework to investigate viscoelastic metadamping in polymeric composite metamaterial beams, in both the infinite and finite configurations. Various phenomena exhibited in the damping-frequency band structures are herein discussed by exploring independent configurations of 1D composites comprising two viscoelastic resonators. For the basis of our comparisons, consider (1) a viscoelastic composite beam hosting two undamped resonators, (2) an undamped composite beam hosting two viscoelastic resonators, and (3) a viscoelastic composite beam hosting two viscoelastic resonators. Later on, we also shed light onto the mechanics of dispersion in high damping regimes as well as the effects of different boundary conditions on the response of the damped beams.

To this end, consider a plain [0 deg] graphite-epoxy composite beam having the same geometric and material parameters given in Sec. 2.3. The lower surface of the beam is attached to a series of 12 equidistant 2-DOF (double) viscoelastic resonators that are *a* = 30 mm apart. Each double resonator consists of two lumped masses that are connected in series to the host beam via viscoelastic springs, as shown in Fig. 1(c). The two discrete masses are selected, such that *m*_{1} + *m*_{2} = *m*_{v}, with *m*_{1} = 0.35*m*_{b} and *m*_{2} = 0.05*m*_{b}, thus ensuring an equal total mass to that of the phononic beam inspected earlier. The first and second viscoelastic resonators are characterized by elastic stiffnesses (*k*_{1},*k*_{2}), viscous damping coefficients (*c*_{1},*c*_{2}), and relaxation parameters (*μ*_{1},*μ*_{2}), respectively. Moreover, the two resonators are tuned such that their natural frequencies $f1=(1/2\pi )k1/m1$ and $f2=(1/2\pi )k2/m2$ are equal with *f*_{1} = *f*_{2} = 1160 Hz.

### 3.1 Influence of Structural Viscoelasticity.

The first configuration consists of a viscoelastic composite beam hosting two undamped resonators (*c*_{1} = *c*_{2} = 0). Consistent with the plain and phononic beams in Sec. 2.3, the same proportionality constant (*α* = 3.58 × 10^{−6}) is used for the hosting beam. Figures 4(a) and 4(b) depict the frequency and damping ratio band structures for a representative unit cell with *μ*_{b} = 3 × 10^{8}. More importantly, these two figures are complemented with close-up insets, revealing band hybridizations that are associated with local resonance bandgaps spanning two different frequency zones (959.8–1072 Hz and 1398–1484 Hz), which lie below and above the resonator’s tuning frequency. Figure 4(c), on the other hand, displays the damping-frequency band structure using two relaxation parameters: *μ*_{b} = 3 × 10^{4} and *μ*_{b} = 3 × 10^{8}, complemented with the modal damping ratios of the finite beam corresponding to the first *μ*_{b} value. The close-up inset provided in Fig. 4(c) reveals that this particular configuration is responsible for producing substantially small damping ratios around the edges of both bandgaps, irrespective of the viscoelasticity level in the beam. Such unique behavior was absent from the phononic beam (see Fig. 3(c)) and is attributed to the local resonance properties of the metamaterial beam. Another critical aspect throughout this effort is how Bloch damping ratios corresponding to the second band (henceforth plotted in green to distinguish it from other bands) change as a function of frequency. To fully understand this, an in-depth study of the viscoelastic effects emerging from the local resonators is required as will be outlined next. Finally, it is worthwhile to mention that for the current configuration, the second band follows a concave-down path.

### 3.2 Influence of Resonator Viscoelasticity.

Second, we consider a configuration where an undamped composite beam (*α* = 0) hosts two viscoelastic resonators. Such configuration can further be subdivided into (1) a metamaterial with a viscoelastic first resonator and an undamped second resonator and (2) a metamaterial where the first resonator is undamped, while the second is viscoelastic. The motivation behind this classification will become apparent shortly. Starting with the first metamaterial, Figs. 5(a) and 5(b) show the frequency and damping ratio band structures when the first resonator is viscoelastic with *ζ*_{1} = 0.003 and *μ*_{1} = 10^{6}, while the second resonator is undamped (*c*_{2} = 0). The local resonance bandgaps induced in this case are found to be identical to the ones corresponding to the first configuration in Sec. 3.1. Figure 5(c) depicts the damping behavior of the metamaterial in the 0–4 kHz frequency range. It shows that damping in the first resonator is mainly responsible for amplifying the Bloch damping ratios around both bandgap edges. More interestingly, the second band separating the two bandgaps now exhibits a concave-up shape. The previous damping features are validated by plotting modal damping ratios of the finite beam in the same figure. Furthermore, unique viscoelastic characteristics of the metamaterial under consideration can be interpreted graphically by computing damping-frequency band structures using different relaxation parameters. For example, Figs. 6(a) and 6(b) correspond to *μ*_{1} = 3 × 10^{3} and *μ*_{1} = 9 × 10^{3}, respectively. A careful inspection of both figures reveals a change in the damping behavior. For *μ*_{1} = 3 × 10^{3}, the maximum attainable Bloch and modal damping ratios of the acoustic band (blue) are larger than their optical band counterparts (red). As such, the second band starts with higher damping ratios than those which appear at the band end. As the relaxation parameter increases to *μ*_{1} = 9 × 10^{3}, this behavior is reversed and the maximum attainable Bloch and modal damping ratios of the acoustic band become smaller than their optical band counterparts, with the second band starting with lower damping ratios than those at its end. In Fig. 6(c), the damping-frequency band structure is computed using several relaxation parameters. It shows that at a certain threshold of the viscoelasticity level in the first resonator, the damping behavior of the resultant metamaterial switches as discussed earlier. Nonetheless, the second band maintains the concave-up shape for all the relaxation parameters used to construct Fig. 6(c).

The mechanism of energy dissipation due to the second resonator is radically different. Figures 7(a) and 7(b) show the frequency and damping ratio band structures, respectively, for a metamaterial beam with an undamped first resonator (*c*_{1} = 0) and a viscoelastic second resonator with *ζ*_{2} = 0.003 and *μ*_{2} = 10^{6}. While some features corresponding to the previous metamaterial (with *μ*_{1} = 10^{6} and *c*_{2} = 0) remain unchanged, such as for example the frequency regions spanned by both bandgaps (depicted in Fig. 7(a)) and the general shape of the acoustic and optical bands in Fig. 7(c), a unique and distinguishable feature of the current metamaterial can be clearly seen from the shape of the second band in Fig. 7(c), which reveals a concave-down formation in contrast to Fig. 5(c). To explore the evolution of the metamaterial’s behavior as a function of *μ*_{2}, damping-frequency band structures are computed in Figs. 8(a) and 8(b) using *μ*_{2} = 10^{3} and *μ*_{2} = 5 × 10^{4}, respectively. Similar observations to the ones made in Figs. 6(a) and 6(b) can be made here as well. In brief, the principal difference between Figs. 8(a) and 8(b) is a dissipative interchanging behavior. Finally, Fig. 8(c) reveals that as *μ*_{2} increases gradually, the resultant metamaterial switches its damping behavior at some point. However, once again, the second bands exhibit a concave-down shape for all the selected values of *μ*_{2}.

Designing such structures with multiple bandgaps can be serviceable in many applications that require vibroacoustic attenuation across broadband frequency ranges or even narrow yet distinct ones. However, a major hurdle is the existence of a number of unavoidable vibration modes between every two consecutive bandgaps. Consequently, a natural question that arises from this discussion is whether it is possible to mitigate such vibration modes without compromising on the efficiency of the induced bandgaps. Previous efforts have already shown that excessive amounts of damping can very well cause bandgap curtailment in locally resonant metamaterials [23]. In particular, incorporating damping in the local resonators without sufficient knowledge can cause severe damage to the metamaterial’s dynamic performance. As such, it is imperative to be conscious of dissipative losses in the local resonators when designing multibandgap metamaterials. To this end, we assemble Figs. 9(a)–9(c), which display the band structures in Fig. 7 superimposed on their corresponding counterparts in Fig. 5. They represent a fair comparison between the metamaterial’s behavior when dissipative effects arise due to the first resonator (alone) or the second resonator (alone), ensuring equal damping amounts in both resonators.

Two main observations can be drawn from Fig. 9(c): (1) while both resonators affect the acoustic and optical bands in the same fashion, amplifying their associated Bloch damping ratios around the bandgap edges, the first resonator is more efficient than the second resonator in attaining such amplification. Moreover, it is clear that damping in the first resonator impacts a wider range of frequencies than the second resonator does. (2) Damping in the second resonator is more authoritative than damping in the first resonator for achieving vibration attenuation between induced bandgaps. The validity of these observations, obtained from unit cell dispersion analyses, can be readily checked by plotting frequency response functions (FRFs) of the corresponding finite beams as shown in Fig. 9(d). Forces with magnitudes of 1 N are applied at 30 mm from the fixed ends of both beams, and the corresponding transverse displacements are calculated at the beam tips. The three close-up insets provided in Fig. 9(d) clearly confirm the main observations stated earlier. However, the effects of resonator damping will be revisited when the time domain analysis is conducted in Sec. 3.3.

### 3.3 Influence of Combined Viscoelasticity.

A more realistic configuration is one which accounts for simultaneous viscoelastic dissipation from the composite beam and the local resonators. For brevity, only damping-frequency band structures will be used to study this configuration as shown in Fig. 10. All displayed subfigures are obtained for a viscoelastic composite beam with *μ*_{b} = 10^{6}. However, each subfigure is associated with a unique resonator damping arrangement. For the first subfigure (Fig. 10(a)), the first resonator is viscoelastic (*ζ*_{1} = 0.1), while the second is kept undamped (*c*_{2} = 0). A number of band structures are computed for relaxation parameters ranging between *μ*_{1} = 2 × 10^{2} and *μ*_{1} = 10^{4}. The results reveal that not all the corresponding second bands have the same shape. For relatively small relaxation parameters (e.g., *μ*_{1} = 2 × 10^{2}), resulting second bands exhibit concave-down formations. However, as *μ*_{1} gradually increases, the shape of the second band changes in return, switching to a concave-up formation at some transition point and maintaining the new shape beyond that point. Changing the damping behavior in this manner implies that the second band shape is controlled by the relative damping of the structure and viscoelastic resonators. In Fig. 10(b), the first resonator is undamped (*c*_{1} = 0), whereas the second is viscoelastic (*ζ*_{2} = 0.1) with relaxation parameters ranging between *μ*_{2} = 2 × 10^{2} and *μ*_{2} = 8 × 10^{3}. Unlike the previous resonator damping arrangement, this one is associated with concave-down second bands for all the considered *μ*_{2} values, which confirms the authority of the second resonator on the second band’s damping behavior, irrespective of the viscoelasticity level in the hosting beam. Finally, when both resonators are viscoelastic with *ζ*_{1} = *ζ*_{2} = 0.1, Fig. 10(c) displays the corresponding band structure for relaxation parameters ranging between *μ*_{1} = *μ*_{2} = 2 × 10^{2} and *μ*_{1} = *μ*_{2} = 7 × 10^{3}. In this case, the combined mechanism of dissipation from all possible damping sources yields the second band formations shown in the figure.

In summary, although the mechanics of dispersion in multibandgap metamaterials are fundamentally controlled by the relative damping of the hosting structure and the resonant inclusions, each resonator damping arrangement is shown to be associated with its own damping characteristics as discussed earlier. For a given finite metamaterial beam, it is interesting to see how different arrangements impact the temporal response in the context of a time-transient analysis. To execute this task, let us consider two finite cantilevered metamaterials that are both assembled from the same viscoelastic beam (*α* = 1.19 × 10^{−5} and *μ*_{b} = 10^{6}). However, it is assumed that the first beam is integrated with a viscoelastic first resonator (*ζ*_{1} = 0.015 and *μ*_{1} = 10^{6}) and an undamped second resonator (*c*_{2} = 0), whereas the second beam is integrated with an undamped first resonator (*c*_{1} = 0) and a viscoelastic second resonator (*ζ*_{2} = 0.015 and *μ*_{2} = 10^{6}). Figure 11(c) shows a comparison between the temporal responses of the two beams when equal forces are applied at 30 mm from their fixed ends. Displacements are calculated at the beam tips when excitations of the form $f(t)=Foe\u2212(t\u2212to)2/2\sigma 2cos(2\pi fct)$ are applied such that the different parameters in the equation are given as follows: *F*_{o} = 500 N, *t*_{o} = 15 ms, *σ* = 1.25 ms, and *f*_{c} = 1265 Hz. The parameters are carefully selected to represent a Gaussian pulse (Fig. 11(a)), which mainly impacts the second band with a spectral content that is centered around *f*_{c} as depicted from the corresponding Fourier spectrum in Fig. 11(b). The comparison shown in Fig. 11(c) validates the findings of Sec. 3.2.

### 3.4 Influence of High Viscoelastic Damping.

As discussed earlier, refraining from using excessive damping in locally resonant metamaterials is often essential to avoid performance deterioration. However, such undesirable performance is yet to be investigated in the context of unit cell predictions. To better understand the mechanics of dispersion within high damping regimes, different energy dissipation sources will be analyzed individually then in combined form. First, only structural damping is increased (*α* = 1.59 × 10^{−4} and *μ*_{b} = 3 × 10^{8}), while both resonators are kept undamped. Figure 12(a) shows that the second bandgap disappears completely although operating in the oscillatory frequency range. The first bandgap is not terminated but is curtailed to a certain extent (it would completely disappear if structural damping is further increased—not shown here for brevity). Moreover, Bloch damping ratios in the third band are shown to increase until reaching the nonoscillatory zone (Fig. 12(b)). Bandgap termination in this particular case (with high structural damping) manifests itself in the damping-frequency band structure in intriguing fashion that can be clearly noticed from the close-up inset in Fig. 12(c). It can be seen that upon reaching a certain point, the second and third bands drop abruptly, heading down in the same direction before splitting in opposite directions at some lower damping point. Such behavior can be realized by superimposing modal damping ratios of the corresponding finite beam on the same figure. The obtained results reveal that as long as the resonators are undamped, the metamaterial beam will still produce frequency regions associated with strong oscillations, regardless of the existence/absence of a frequency bandgap and irrespective of the damping level in the hosting beam. Conversely, when the damping level in the resonators is increased while the host beam is kept undamped, the resulting dispersion characteristics substantially change. Consider a case where only the first resonator is highly damped (*ζ*_{1} = 0.6 and *μ*_{1} = 10^{6}), while the second resonator is undamped. Figures 13(a) and 13(b) show that all three bands intersect with each other, terminating both bandgaps even before reaching the nonoscillatory zone. A notable feature of Fig. 13(c) is the mechanism of the second bandgap termination, which shows that on reaching a certain point, the second and third bands rise abruptly, heading up in the same direction before splitting in opposite directions at some higher damping point. This behavior can be realized by superimposing modal damping ratios of the corresponding finite beam on the same figure. Finally, the combined effect of high structural and resonator damping is shown in Fig. 14 for three different cases: (1) moderate structural damping combined with high first resonator damping as shown in Fig. 14(a), (2) high structural damping combined with high first resonator damping as shown in Fig. 14(b), and (3) light structural damping combined with high damping in both resonators as shown in Fig. 14(c). All three cases demonstrate the complete extinction of frequency bandgaps as a result of cumulative damping. Dispersion characteristics similar to the ones deduced from Figs. 12 and 13 are also observed in Fig. 14. It should be noted that these cases represent only a few examples of an infinite number of possible damping combinations. Nevertheless, we believe the results presented here describe the most important dispersion features encountered in locally resonant metamaterials.

### 3.5 Influence of Physical Boundary Conditions.

All the infinite medium predictions have thus far been validated using finite beams with only one type of physical boundary conditions, namely, fixed-free (or cantilever). We briefly demonstrate the robustness of such unit cell predictions by considering finite beams with three different physical boundary conditions: (1) hinged-hinged, (2) fixed-hinged, and (3) fixed-fixed. In Fig. 15, the modal damping ratios of each finite beam are computed and superimposed on the damping-frequency band structure obtained from a representative unit cell to produce three subfigures similar to the one provided in Fig. 9(c). Interestingly, vibration modes close to bandgap regions exhibit damping ratios that are only very minimally affected by the type of boundary conditions. These results confirm that the most important metadamping features are predominantly controlled by the presence of the resonant inclusions in multibandgap metamaterials.

## 4 Design for Multiple Bandgaps

To provide general guidelines for the design of multibandgap viscoelastic metamaterials, consider an undamped composite beam (*α* = 0) leveraging triple serially connected local resonators. To understand the sole damping effect from each resonator, Fig. 16 is created. Each of its three subfigures depict the metadamping mechanics when two resonators are kept undamped, while the third is rendered lossy with *ζ* = 0.003 and *μ* = 3 × 10^{8}:

When only the first resonator is damped and the second and third are kept undamped (

*c*_{2}=*c*_{3}= 0), resulting intermediate bands (second and third) produce the concave-up formations shown in Fig. 16(a). The previous implies that this particular resonator damping arrangement produces substantially small Bloch damping ratios at some intermediate points through each of these two intermediate bands.When only the second resonator is damped while the first and third are kept undamped (

*c*_{1}=*c*_{3}= 0), the second band ends with substantially small Bloch damping ratios, while the third band begins with relatively small Bloch damping ratios (*ζ*< 0.5 × 10^{−3}) as shown in Fig. 16(b). Moreover, this resonator damping arrangement is less efficient than the previous one in terms of attenuating vibration modes belonging to the acoustic (blue, leftmost) and optical (red, rightmost) bands.When only the third resonator is damped while the first and second are kept undamped (

*c*_{1}=*c*_{2}= 0), relatively small Bloch damping ratios are observed at the start of the second band (*ζ*≅ 0.7 × 10^{−3}) as well as at the end of the third band (*ζ*≅ 0.5 × 10^{−3}) as shown in Fig. 16(c). Furthermore, this resonator damping arrangement is the least efficient in attenuating vibration modes belonging to the acoustic (blue, leftmost) and optical (red, rightmost) bands.

To achieve broadband vibroacoustic attenuation that extends from the beginning of (and slightly before) the first bandgap to the end of (and slightly after) the third bandgap, it is clear that between the three arrangements discussed earlier, the third is the best choice. However, such damping performance can be further enhanced by utilizing an arrangement that incorporates simultaneous and concurrent damping in both the second and third resonators. To ensure a fair comparison, total damping (*ζ* = 0.003) is equally divided between the second and third resonators such that *ζ*_{2} = *ζ*_{3} = 0.0015 with *μ*_{2} = *μ*_{3} = 3 × 10^{8} to produce the damping-frequency band structure shown in Fig. 17, which is superimposed on its counterpart from Fig. 16(a). Upon careful inspection of this figure, it is clear that the second and third bands produce concave-down formations that exhibit Bloch damping ratios >0.001 without compromising on the efficiency of the induced bandgaps.

## 5 Conclusions

An in-depth analysis of the elastic wave dispersion mechanics was conducted in this article to investigate the viscoelastic damping behavior of infinite as well as finite 1D polymeric composite metamaterials with multiple bandgaps. The applied free wave unit cell model was derived based on well-established theories in linear viscoelasticity and wave propagation in periodic solids. The proposed framework was validated against commonly used benchmark viscous damping models. Following this, the model was deployed to elucidate the viscoelastic metadamping phenomenon in a graphite-epoxy metamaterial beam comprising periodic double and triple resonators. It was concluded that the dispersion mechanics associated with a given metamaterial are highly sensitive to the chosen resonator damping arrangement. Specifically, for a double resonator metamaterial beam, it has been found that while both resonators affect the acoustic and optical bands in the same fashion, amplifying their associated Bloch damping ratios around the bandgap edges, the first resonator is more efficient than the second resonator in attaining such amplification. Moreover, damping in the first resonator was shown to be capable of affecting a wider range of frequencies than the second resonator. Conversely, damping in the second resonator is more authoritative than damping in the first resonator for achieving vibration attenuation between induced bandgaps. Depending on the selected resonator damping arrangement, damping-frequency band structures generate intermediate bands with different shapes and concavities, the understanding of which is critical to the effective design of bandgap structures for applications that require vibroacoustic attenuation across broadband frequency ranges (or even narrow yet distinct ones). Finally, the robustness of the unit cell predictions has been demonstrated. Results show that vibration modes close to bandgap regions exhibit damping ratios that are only minimally affected by the type of boundary conditions, implying that the most important metadamping features are predominantly controlled by the presence of the resonant inclusions in multibandgap metamaterials.

## Acknowledgment

M. Nouh acknowledges the support of this work from the US National Science Foundation through Award no. 1847254 (CAREER).

## Conflicts 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 authors upon reasonable request. The authors attest that all data for this study are included in the paper.