## Abstract

The study of phononic materials and structures is an emerging discipline that lies at the crossroads of vibration and acoustics engineering and condensed matter physics. Broadly speaking, a phononic medium is a material or structural system that usually exhibits some form of periodicity, which can be in the constituent material phases, or the internal geometry, or even the boundary conditions. As such, its overall dynamical characteristics are compactly described by a frequency band structure, in analogy to an electronic band diagram. With roots extended to early studies of periodic systems by Newton and Rayleigh, the field has grown to encompass engineering configurations ranging from trusses and ribbed shells to phononic crystals and metamaterials. While applied research in this area has been abundant in recent years, treatment from a fundamental mechanics perspective, and particularly from the standpoint of dynamical systems, is needed to advance the field in new directions. For example, techniques already developed for the incorporation of damping and nonlinearities have recently been applied to wave propagation in phononic materials and structures. Similarly, numerical and experimental approaches originally developed for the characterization of conventional materials and structures are now being employed toward better understanding and exploitation of phononic systems. This article starts with an overview of historical developments and follows with an in-depth literature and technical review of recent progress in the field with special consideration given to aspects pertaining to the fundamentals of dynamics, vibrations, and acoustics. Finally, an outlook is projected onto the future on the basis of the current trajectories of the field.

## Introduction

A periodic medium is a material or structural system that exhibits some form of spatial periodicity, which can be in the constituent material phases, or the internal geometry, or the boundary conditions. The study of periodic materials and structures has a long history in the field of vibrations and acoustics with origins that can be traced back to Newton's first attempt to describe the propagation of sound in air [1] and Rayleigh's early study of continuous periodic structures [2]. The topic has grown to become of fundamental importance in condensed matter physics due to the role that atomic vibrations (and electronic structure) play in determining the properties of crystals. The notion of a “phonon” has emerged in the context of vibrations in a crystal lattice. Formally defined as a quantum of vibrational energy in an elastic medium—which may be interpreted as a discrete particle-like quantity of sound in a solid—the term has also been associated with classical wave-like vibrations and acoustics, mainly in the context of periodic media. Subsequently, it has grown to become customary to refer to a periodic material, or structure, even if in the size scale of a large engineering system, as a “phononic material,” or a “phononic structure.” As such, there has been an exquisite connection, especially in the past few decades, between phonon physics and dynamics of periodic materials and structures. Perhaps what stands out the most is the concept of a band diagram—a diagram that represents the relationship between frequency (or energy) and wavenumber, along multiple directions. In physics, this diagram represents the backbone of electronic structure theory, credited for forming a basis for the classification of all crystals into metals, semiconductors, and insulators. In mechanics, a band diagram is precisely a representation of the dispersion relation describing the nature of free wave propagation in an elastic (or acoustic) medium.

Looking closely at the study of periodic systems in the past half-century, we find that researchers in vibrations and acoustics, and more broadly from the mechanics community at large, have conducted a considerable amount of work on key theoretical foundations, concepts, and analysis techniques that are relevant to periodic systems in other nonmechanics disciplines. Arguably the two most motivating applications in mechanics, going back to the 1950s, and extending through the 1990s and beyond, have been composite materials (which conveniently have been modeled as periodic materials) [3,4] and aircraft structures (which naturally exhibit some degree of periodicity emanating from the presence of ribs introduced primarily for strengthening and stiffening purposes) [5,6]. Other areas related to mechanical and civil engineering include multiblade turbines [7–9], impact resistant foam/cellular materials [10–12], periodic foundations for buildings [13,14], and multistory buildings and multispan bridges [15].

In recent years, starting in the early 1990s, the field has experienced a resurgence with the introduction of phononic crystals [16–21], and, roughly a decade later, with the study of acoustic metamaterials [22], which may be also broadly considered as examples of phononic materials. A phononic crystal is a composite or nonuniform material consisting of one, two, or more material phases (solid and/or fluid) arranged periodically in space. This is not much different from periodic composite materials studied earlier in the engineering literature, with the only difference that more rigor is applied in the treatment of symmetry by borrowing concepts from crystallography [23]. The already existing connection between periodic materials and phonon physics has received a substantial upgrade with the introduction of phononic crystals. An acoustic (or elastic) metamaterial, on the other hand, is also (usually) a periodic material except it has the added feature of exhibiting local resonance, or possibly other traits, that can give rise to “unusual” dynamical behavior, such as negative refraction and negative group velocity. Metamaterials also have a foundation in electromagnetics, and the disciplines of acoustics and elastodynamics have been steering a parallel path in the development of this highly promising class of materials.

While applied research in phononic materials has been abundant in recent years, inspection from a fundamental mechanics perspective, and particularly from the standpoint of dynamical systems, is needed to advance the field in new directions. For example, techniques already developed for the treatment of damping and nonlinearities in structural dynamics provide a rich resource for application to the problem of wave propagation in periodic materials and structures. Similarly, numerical and experimental approaches originally developed for the characterization of conventional materials and structures provide indispensable tools that may be employed toward better understanding and exploitation of phononic systems.

The goal of this review article is to bring the topic of phononic materials to the attention of the Journal of Vibration and Acoustics and the Applied Mechanics Reviews communities, and to encourage interest in the field among researchers and practitioners from across all related disciplines. To this end, the article begins with an overview of historical developments underlying the early years of research in periodic materials and structures, and follows with an in-depth literature review of progress that came about with the advent of phononic crystals and acoustic/elastic metamaterials with special consideration given to aspects pertaining to dynamics, vibrations, and acoustics. A detailed technical review is then presented on fundamental concepts in phononic materials theory, numerical techniques, and the treatment of damping and nonlinearity, as well as experimental research encompassing the field. Finally, a brief outlook is projected toward the future on the basis of the progress achieved to date.

### Historical Origins.

According to Ref. [24], the first work on a periodic one-dimensional (1D) lattice is the result of Newton's attempt to derive the formula for the velocity of sound in air. One can speculate that because of lack of knowledge of differential calculus, a natural approach would have consisted in discretizing the continuous media supporting the propagation of sound into a series of lumped masses connected by a lumped spring. This process, which is illustrated in Fig. 1, leads to a periodic system, which is known to have a more complex wave behavior than the original continuous system.

Fig. 1
Fig. 1
Close modal

The same system was later investigated by John Bernoulli and his son Daniel, who, in a series of studies starting in 1727, demonstrated that a system of N masses is characterized by N modes of vibration and associated frequencies, and essentially formulated the principle of superposition. Subsequent studies included the estimation of the velocity of wave propagation along one axis of a cubic lattice structure as a function of wavelength. These efforts have led to the first instance of observation of wave dispersion from the estimation of the phase velocity's dependency on frequency. These results were later expanded and explained in 1881 by Lord Kelvin [25] among others, and the reader is referred to the classic books by Brillouin [24] and Kittel and Kittel [23] for a detailed historical perspective on the topic.

### Periodic Materials and Structures in the 20th Century

#### Composite Materials.

Periodic layered composites are the motivation for one of the first studies, as in Ref. [26], where a 1D periodic arrangement of different materials is investigated. The study of this class of assemblies found much interest in the 1960s and 1970s with the advent of composite materials, and the necessity to assess their integrity through nondestructive evaluation (NDE) methods based on the propagation of elastic waves. This has motivated studies of layered composites and the development of equivalent, or continuum models, as in the case of Refs. [3,27]. A continuum theory was developed therein to describe the dynamic behavior of a laminated composite by expressing the displacements of the reinforcing and the matrix layers as two-term expansions of the layers' midplanes. The accuracy of the theory was estimated by evaluating the dispersion curves for waves with polarizations normal and parallel to the layering and comparing them with exact curves. Variational formulations for the analysis of harmonic wave propagation in periodic elastic composites were presented by Nemat-Nasser in Ref. [4], which is the first paper of a series presenting studies on layered configurations of increasing complexities [28–32]. An early review of variational methods for the analysis of periodic composites can be found in Ref. [33], while continuum theories have also been proposed by Nayfeh and coworkers in the same period [34]. Their study of wave propagation in layered composites, again mostly motivated by NDE applications, has continued through the years, and notable contributions can be found in Refs. [35,36] along with studies on guided waves in layered composites as in the case of Refs. [37–39]. In several studies, the transfer matrix approach and other, similar matrix-based calculations are applied to investigate dispersion properties and the propagation of harmonic waves through the thickness of composites of general in-plane anisotropy. Mixture theories for the estimation of equivalent stiffness and mass properties of layered structures have also been pursued by other authors; see, for example, Refs. [40,41]. The ability of equivalent, mixture, or continuous theories to predict the dispersion properties of heterogeneous media is still the object of intense investigation, and has been addressed in some of the multiscale studies performed in recent years. Examples are the studies by Hulbert and coworkers [42–44], and by Fish and others [45–47] who have considered the 1D periodic waveguide as a test case for the initial assessment of the developed approximation techniques based on multiple scale expansions. Periodic configurations of highly mismatched materials have formed the basis for most of the periodic systems considered by the community, and are still a common trait in the design of “phononic crystals.” One of the early investigations in this particular area is found in Ref. [18], where two-dimensional (2D) periodic materials and their complete bandgap properties are discussed as a parallel to what was previously observed in “photonic crystals” (see Sec. 1.3).

#### Periodic Structures.

In the area of structural engineering, Cremer and Leilich [48] were among the first to investigate harmonic flexural wave motion along a 1D periodic beam. The periodicity was introduced by either simple supports or by point masses located at regular intervals. The periodic beam is an example of a monocoupled periodic system whereby each unit cell communicates with its neighbors through a single degree of freedom. Its wave characteristics at any frequency are therefore described by a single wave mode and by a single pair of equal and opposite propagation constants. Other examples of earlier studies are reported in Ref. [49] and in the works of Heckl [50,51] and Ungar [52]. Therein, investigations were conducted on beam-plate systems, 1D beam structures with various periodic features and 2D beam grillages. Most of these systems and their propagation characteristics can be conveniently investigated through the formulation of a transfer matrix. Its repeated application through the number of cells composing a finite periodic system leads to the estimation of natural frequencies and mode shapes. In addition, the analysis of a single unit allows the estimation of the dispersion relations according to a procedure briefly summarized later in this paper (see Sec. 2.2). Among the pioneering efforts on the formulation of the transfer matrix, the work of Lin and McDaniel [53] is worthy of note, along with the work of Faulkner and Hong [54], which provides examples that can be useful to researchers who want to begin investigations in this area. Of interest is also the work of McDaniel [55,56], Thomas [57], and Williams [58] on rotationally periodic structures, which is an area that appears to have developed on an independent path, mostly as a result of its strong practical relevance and associated specific technological problems. Prime examples of rotationally periodic structures are bladed disks in turbomachinery, where several characteristic dynamic properties of periodic structures have been observed and linked to large responses and localization of vibrations, and have been identified as causes of failure [59,60].

Many noted contributions in this area are attributed to Mead and his coworkers at the University of Southampton, which has been a hotbed of activities starting from the 1960s. An excellent review of their work has been produced by Mead himself in Ref. [61], and will not be repeated here for the sake of brevity. Among the Southampton studies, we would like to mention the work of Gupta [62,63], which relates the modal content of beams, stiffened panels, and ribbed plates to their dispersion properties. The investigation of modal properties of monocoupled and multicoupled periodic systems can also be found in Refs. [64,65]. Reference [64], in particular, relates the bounding frequencies of the propagation zones of monocoupled systems to the natural frequencies of the individual unit cells and of a finite assembly. The work is extended to multicoupled systems in Ref. [65], where the effect of damping is also considered. The forced response of periodic systems is studied in Refs. [66,67], where a response function is found for an infinite, uniform, 1D structure subjected to an array of harmonic forces or moments. The case of random convecting loading is addressed in Refs. [68,69]. The application of numerical techniques for the prediction of dispersion and overall dynamic response is presented in Ref. [70], where the finite element (FE) implementation of wave periodicity conditions as defined by Bloch's theorem is illustrated and applied to simply supported beams and to complex built-up structures. The formulation of a real-valued eigenvalue problem for dispersion estimation, which is suitable for its solution through built-in modal solvers in commercial FE packages, was later presented by Thomas [57] and subsequently by Åberg and Gudmundson [71]. The application of the FE method for the analysis of periodic structures is an area of continued interest, as indicated, for example, by the work of Mace and coworkers [72,73]. In particular, the periodicity of the structure, or the consideration of a uniform structure as the limiting case of a periodic assembly, can lead to significant advantages in terms of computational cost, or to the ability to exploit commercial FE codes for the derivation of the dispersion properties of complex structural components, such as stiffened plates, sandwich panels, or cylindrical shells [73]. These approaches have been successfully extended to the study of 2D structures, as in Ref. [74]. The analysis of 2D periodic structural components, with a focus on their directional properties in addition to the stop and pass behavior, traces back to the early studies of Langley (see, for example, Refs. [75–77]) and continues to draw interest in regards to the influence of topology and geometry on the wave propagation characteristics of structural lattices [78–83]. The presence of forbidden zones of response, the observation of “caustics” [84], and the relation between wave directionality and the equivalent anisotropic elastic properties of the lattice are of relevance to the design of impact absorbing structures, acoustic waveguides, vibration isolation devices, and solids with tailored mechanical properties.

### Phononic Crystals.

As mentioned earlier, a remarkable trait of phononic materials/structures research is the broad interest across disciplines. While research in composite materials and periodic structures has evolved among mechanicians the notion of artificial periodicity has also been appealing to the electromagnetics and photonics communities (see Ref. [85] for a review of historical developments in these areas). A key milestone occurred in 1987 when Yablonovitch [86] and John [87], almost simultaneously, proposed the concept of what was later referred to as a photonic crystal—a multiphase dielectric material ordered with symmetries resembling those of atomic-scale crystalline materials (as dictated by the rules of crystallography). Only a few years later, the same concept emerged for acoustic or elastic waves, which was subsequently called a “sonic crystal” or phononic crystal [16–21]. The public had a luminous introduction to the concept when a large periodic array of cylinders was put on display in Madrid to represent a sound sculpture [88].

#### Configurations of Phononic Crystals.

In 1992, Sigalas and Economou [16] considered in-plane (mixed longitudinal-transverse) elastic wave propagation in a 2D phononic crystal, albeit composed of fictitious materials. In the following year, they repeated their study using gold cylinders embedded in a beryllium host and incorporated out-of-plane shear waves [17]. Kushwaha et al. [18,19], on the other hand, studied out-of-plane waves in a periodic array of nickel alloy cylinders in an aluminum alloy matrix. An analysis of in-plane wave propagation in a 2D phononic crystal consisting of carbon cylinders as inclusions and an epoxy resin as a matrix material was conducted around the same time by Vasseur et al. [20]. The problem of a three-dimensional (3D) phononic crystal consisting of solid sphere inclusions in a solid host was also considered by Economou and Sigalas in 1994 [89], and later by Psarobas and coworkers [90] and Sigalas and Garcia [91]. Meanwhile, there has also been significant interest in sonic crystals consisting of multiple phases of liquids and/or gases. The acoustic band structure of these systems were first reported by Sigalas and Economou [17], followed by other studies, such as the theoretical investigations by Kushwaha and Halevi [92,93] and the experimental investigations by Sánchez-Pérez et al. [94] (audible range) and Montero de Espinosa [95] (ultrasonic range). Regardless of the phase state and choice of materials, two distinct features were identified as critical for the design of phononic crystals: the unit-cell topology and the lattice symmetry. It was shown that for elastic waves a “cermet topology,” in which a stiff/dense material forms the matrix and a compliant/light material is used as inclusions, generally produces the largest bandgaps [96]. More recently, this concept was taken further by virtue of full-fledged unit-cell topology optimization [97,98]. As for lattice symmetry, it was revealed that phononic crystals patterned over a hexagonal lattice typically display larger bandgaps compared with those based on square lattices [99]. The specific choice of the constitutive material phases and the possibility of interweaving more than one lattice also provide means for improving bandgap characteristics [100].

In the late 1990s, the concept of phononic crystals was extended to free surfaces involving 2D semi-infinite domains consisting of a periodic array of inclusions extending indefinitely into the depth in the direction normal to the surface [101–104]. In this configuration, Rayleigh-type Bloch waves propagate along the surface. Among numerous investigations, the case of anisotropic material properties was considered in this context [101,104]. Similar to free surfaces, phononic crystals were also studied and realized in the form of plates of finite thickness [105–117]. The interest in this class of problems has been mostly in the study of Rayleigh–Lamb waves in plates with cylindrical inclusions composed of either a second solid phase [106,107,109] or a fluid, such as air or even simply a vacuum [112–114]. The possibility of using spheres as inclusions in plate-based phononic crystals has been studied as well [108]. Contribution of phononic crystals to silicon MEMS [118] and standard MEMS [119] has emerged from this line of research. Further analysis and physical realization of phononic crystal plates with unit-cell sizes at the micrometer scale created a platform for radio-frequency sensing applications [112–114]. For additional detail, the reader is referred to the review article by Olsson and El-Kady [115], and to other recent reviews on phononic crystal plates, or 2D phononic crystals in general, which can be found in Refs. [116,117].

As research in phononic crystals continued to grow, other configurations were proposed, for example, by Torrent and Sánchez-Dehesa, who considered a radial phononic crystal for which a frequency band diagram was calculated for wave motion in radial coordinates [120]. More subtle symmetries have also been considered, such as quasicrystal [121] and fractal [122,123] phononic crystals—in analogy to similar studies concerned with photonic crystals. Furthermore, the field also witnessed advances in the type of theoretical investigations carried out. For example, Laude et al. [124] and Romero-García et al. [125] devised techniques to compute evanescent Bloch wave modes associated with phononic crystals alongside the usual propagating modes. The effects of periodicity truncation rendering phononic crystals finite in their spatial extent have also been examined, e.g., see [126–129], which builds on earlier studies on the topic by Gupta et al. [62,63].

#### Applications.

The realization of phononic crystals in different forms, e.g., bulk, semi-infinite surfaces, plates, beams, has opened up a range of opportunities for realistic applications in engineering and applied physics. The concept of a bandgap naturally lends itself to applications involving vibration and acoustic filtering and control. For example, vibration isolation [130,131] or minimization [128,132], as well as noise isolation and control [94,133], may be realized with profound effectiveness. Similarly, nondestructive testing of materials for defects may be conducted by utilizing phononic crystal phase-surface characteristics for the emission of probing “wave beams” [80].

One of the remarkable functions that was realized using phononic crystal is waveguiding and localization by means of defect-engineering [134–139]. This has been achieved by introducing line or point “defects” into an otherwise perfect phononic crystal. The defect is incorporated by simply removing one or more inclusions (or other features that the characteristic unit cell possesses) in a manner as to introduce an anomaly in the form of a given pattern. In the defected region, the system will enable mode localization within the phononic crystal's bandgap frequencies, and this in turn allows for waveguiding or localization in a cavity. The latter phenomenon provides the basis for signal sensing in the microscale phononic crystals mentioned above [112–114]. These systems hold promise for radio-frequency communications due to the practicality of direct integration over CMOS circuitry with low insertion losses and their ability to provide frequency selection with high quality factors. In addition, the principle of defect-based elastic wave confinement has been utilized toward the realization of wave demultiplexers—passive devices that enable waves at certain frequencies to travel in different directions [140–142]. Further development along this track has started to produce early realizations of digital logic circuits based on phononic crystals. For example, phase differences associated with lack of collinearity between the wavevector and the group velocity (as revealed by equifrequency contour diagrams) have been used toward the creation of basic logic gates [143,144].

Utilization as sensors for the determination of liquid properties has been yet another application for phononic crystals [145,146]. The underlying concept is based on introducing a line defect that acts as a slit cavity, which in itself is also part of the fluidic system. The cavity mode causes a distinct transmission peak within the bandgap of the phononic crystal where the frequency of maximum transmission depends on the speed of sound of the confined liquid. In other work involving microfluidics, Cooper et al. [147] have demonstrated a simple interface between a piezoelectric surface acoustic wave device and a disposable microfluidic chip, patterned with a phononic crystal structure. They have illustrated that microfluidic manipulation, such as centrifugation of blood, may be performed on a disposable phononic chip. Employment of phononic crystals for mass sensing has also been demonstrated, e.g., Nardi et al. used a hypersonic surface phononic crystal to realize ultrahigh sensitivity mass sensing at operational acoustic frequencies that exceed 100 GHz [148].

Another active line of applied research is concerned with wave collimation and refraction. The interest here has been primarily in acoustics, although the ideas are also applicable to elastic waves. The function of sound collimation has been studied by Chen et al. [149], Espinosa et al. [150], Christensen et al. [151], and Shi et al. [152]. Refraction and focusing of sound waves have been a major topic due to the potential significant impact on the fields of subsurface imaging and NDE. Phononic crystals can be used to induce two types of refractions: positive and negative. Positive refraction is carried out by shaping the phononic crystal boundaries in order to refract an incoming wave [153–155]. The realization of negative refraction on the other hand is more subtle as it is based on tuning the phononic crystal to exhibit opposite phase and group velocity signs at the frequency of interest [156–163]. Use of phononic crystals for focusing has also been realized using nonlinear solitary waves generated in an array of precompressed metallic beads [164].

Another intriguing application of phononic crystals is wave rectification, whereby a device is designed to allow acoustic or elastic waves to travel along one direction but not in the opposite direction. This principle promises to be of great benefit for vibration mitigation but also could be used as a building block in acoustic/elastic logic gates in analogy to electronics. While not yet perfectly implemented, some promising approaches have been proposed. For example, Liang and coworkers utilized 1D layered phononic crystals composed of two segments, a linear segment to filter an incoming wave and a nonlinear segment to enable energy to transfer to a different mode at a slightly different frequency [165,166]. Nonlinearity was used once again, now in the context of bifurcation phenomenon within a linear chain of precompressed beads, to generate a rectification effect [167]. Two-dimensional linear phononic crystals were employed as well to realize an acoustic diode, but with a mechanism based on diffraction and modal conversion at bandgap edges [168]. A major challenge that is yet to be tackled is to achieve rectification in a perfectly asymmetric fashion, without alteration of frequencies and along a straight line and a broad frequency band.

The coupling of acoustic or elastic waves with other types of waves, such as electromagnetic signals provides yet another avenue for promising phononic crystal applications. For example, a crystal may be designed to exhibit simultaneous phononic and photonic bandgaps coinciding at the same range of wavevectors [169–174]. This can be used to induce an interaction between the two types of waves in a defect-type cavity [175]. This coupling could in turn be exploited for signal manipulation or as a mechanism to slow down light. The simultaneous phononic/photonic bandgap concept is now often linked with the term “phoXonic crystal,”2 and one of the first investigations of the concept can be found in Ref. [169]. The concept of wave coupling has been utilized in the context of acousto-optics as well [176]. A recent review article discusses developments in this area as well as other applications of phononic crystals and related concepts [177].

All the work reviewed thus far on phononic crystals (and other types of periodic media) has been based on the concept of passive control of wave propagation. While very practical in its own right in many applications, a natural extension is the development of techniques for active control or “tunability.” Such a prospect certainly holds promise as it implies an “on the fly” capability to change the location and size of bandgaps, as well as the shapes of the dispersion curves themselves. A simple approach for realizing active control of phononic crystals involves anchoring the unit-cell inclusions on axes and thus allowing them to be rotated on command [178]. As such, the orientation of the inclusions with respect to the incoming waves may be tuned thus altering the band structure during operation. A variety of innovative approaches have been proposed to achieve the same objective but by other means. Piezoelectricity, for example, has been utilized with noticeable effectiveness [179–182]. Other approaches involve the use of light-induced substrate potential in colloidal phononic crystals [183] or the use of electrorheological materials in which tuning is done via an external electric field [184]. Magnetism, in its part, has also been considered by several researchers as a platform for realizing phononic crystal tunability through alteration of material elastic properties. For example, Wang et al. [185,186] employed piezomagnetism and examined magnetoelectroelastic phononic crystals, and Robillard et al. [187] considered magnetoelasticity. Temperature control has provided yet another approach for tunability. This was explored as early as 2000 by Ruzzene and Baz [188] in the context of shape memory alloys whose properties are temperature dependent. Another phononic crystal study involving temperature investigated its effect on the bandgaps of surface and bulk acoustic waves [189]. Similarly, Jim et al. [190] considered the possibility of thermal tuning of a phononic crystal composed of a ferroelectric ceramic and epoxy. Finally, we also report on efforts that utilize finite elastic deformation [191] and lattice-level mechanical instabilities [192], both being approaches that take advantage of large elastic deformation.

While the above applications appear at a variety of scales, they all may be practically viewed as macroscopic in nature (e.g., may be modeled using continuum theory). The concept of a phononic crystal, however, can in principle be also realized at the atomic scale whenever the conditions for coherent phonon transport exist. Under such conditions, the dispersion characteristics of both the constitutive material lattice (defined by a primitive cell) and the phononic crystal lattice (now defined by a nanoscale supercell) contribute to the physics of the problem, which mostly pertain to transport properties, such as the thermal conductivity. This viewpoint in itself has opened up an entire track of research that overlaps with the emerging discipline of nanoscale heat transfer [193,194]. The reader is referred to Refs. [195–198] and [199–201] for, respectively, theoretical and experimental investigations of “nanophononic crystals.” More recently, a new analogy between the acoustic/elastic regime and the thermal transport regime has been introduced through the proposal of a “nanophononic metamaterial” whereby thermal conductivity reduction is achieved by local resonance [202]. These two concepts of nanoscale phononic materials provide promising opportunities for advancing a broad range of applications in heat management, such as those discussed in Ref. [203].

### Acoustic/Elastic Metamaterials.

As mentioned earlier, the scope of interest in periodic materials received yet another surge in the early 2000s, and this has not been any less significant than the one brought about by phononic crystals in the early 1990s. The new phononic material that emerged has been denoted acoustic, or elastic, metamaterial, in analogy to an electromagnetic metamaterial. The prefix “meta” in this context is obtained from the Greek language to indicate the notion of “after” or “beyond.” This is linked to the understanding that metamaterials exhibit properties that are generally viewed to go beyond what we expect to find in naturally occurring, or conventional, materials. Examples of such unusual properties include the generation of subwavelength bandgaps and negative refraction. A key feature of acoustic/elastic metamaterials is the presence of local resonance, or possibly other traits that may bring about similar effects. While periodicity is not necessary for an acoustic/elastic metamaterial, its presence enables intrinsic, unit-cell-based, description of dynamical properties and wave propagation characteristics as well as order in the introduction of the locally resonant elements.

#### Configurations of Locally Resonant Metamaterials.

A seminal paper that has introduced a conceptual realization of an acoustic metamaterial was published in 2000 by Liu et al. [22]. In this paper, a 3D array of lead spheres coated with a 2.5-mm layer of silicone rubber was stacked in a simple cubic arrangement within an epoxy matrix. The lattice constant of the periodic material system was 1.55 cm. Upon excitation with acoustic waves, it was revealed that the medium exhibited a bandgap at a wavelength well below the regime corresponding to bandgap generation based on spatial periodicity of the impedance mismatch areas, which leads to “Bragg scattering.” In addition, the localized resonant structure was shown to cause the material to behave as if its effective elastic constants were negative at certain frequency ranges (precisely at frequencies around the bandgaps associated with the local resonances). Interestingly, the form of a locally resonant bandgap has appeared earlier in the Sigalas and Economou [16] paper although not for the purpose of exploring a new bandgap mechanism and the results were based on fictitious materials.

The use of local resonance as a mechanism for the generation of bandgaps and the achievement of other intriguing subwavelength properties has effectively opened up a new major thrust in phononic materials research. After the initial proposal of a three-phase composite material [22,204], numerous other configurations for locally resonant periodic materials have been proposed. These include the use of binary materials with the “inclusion” material being very soft [205,206], hollow cylinders or spheres [207], split rings or spheres [208,209], beams or plates with suspended masses [210–212], plates or surfaces with pillars [213–216], and structures with inertial amplifiers [217].

The notion of a locally resonant bandgap has been a source of studies concerned with analysis of the nature of the bandgap and how it differs from the more familiar Bragg scattering bandgap. Goffaux et al. [218], for example, have shown that the transmission spectrum corresponding to a locally resonant bandgap displays fano-like features. Xiao and coworkers [219,220] provided closed-form analytical expressions for the prediction of bandgap edges, and Liu and Hussein [221] provided a criterion for the transition between Bragg scattering and local resonance as the resonator properties are varied for beam structures with a periodic array of suspended masses.

#### Effective Properties.

Like their electromagnetic counterpart, acoustic/elastic metamaterials enjoy the allure of exhibiting unique dynamical behavior at long wavelengths, which is attributed to the presence of subwavelength bandgaps. This trait has made this type of materials an excellent candidate for the realization of desired homogenized properties at the long-wave end of the spectrum. Consequently, there has been a growing interest in deriving effective properties of locally resonant phononic materials.

Using the coherent potential approximation method [222], Li and Chan [223] derived the effective elastic modulus and effective density of a periodic composite material consisting of soft rubber spheres suspended in water (a choice motivated by the rubber's much slower speed of sound compared with that of water). They showed that at certain frequency ranges near the local resonance frequency, the effective bulk modulus and density become negative. Furthermore, their results have shown that there is a narrow frequency window at which both properties are simultaneously negative. The authors discussed how double negativity relates to rather exotic properties, such as negative refraction. They also correlated the observed phenomenon with the type of resonance modes in the local unit-cell structure. Specifically, it was explained that a negative elastic modulus is associated with a monopolar resonance, and a negative density is associated with a dipolar resonance. Upon reflection on the intriguing phenomenon of negative effective properties, one may interpret these as dynamical conditions whereby the restoring and inertial forces act in the direction of motion (rather than opposite to it) for the cases of effective elastic modulus and effective density, respectively.

Other studies followed with the aim of calculation, analysis, and/or utilization of negative effective properties, demonstrating negative effective elastic modulus [224–227], negative density (or mass) [228–232], or both [233–236]. A natural extension to these is the possibility of anisotropy in the effective mass or density [155,237–241], which is yet another intriguing concept in its own right. While most effective property calculation techniques are based on unit-cell analysis, it is also possible to compute these properties from the reflection and transmission coefficients that may be obtained by simulations or experiments [242]. Furthermore, while most studies have focused on actual, realistic models of locally resonant acoustic or elastic metamaterials, much can be learned by exploring the problem in the context of a simple “mass-in-mass” lumped parameter model as done by Huang et al. [243].

The topic of effective properties of acoustic/elastic metamaterials cannot be viewed in isolation from the classical theory of homogenization (e.g., see Ref. [244] and subsequent high-order homogenization formulations for solids [42–47], and homogenization of rigid bodies in a fluid [245]). However, in order to adequately obtain frequency-dependent homogenized properties, it is necessary to fully generalize the theory of homogenization to the dynamic regime and to do so in a manner that allows for an accurate regeneration of the dispersion relation at any wavenumber and frequency. This problem, commonly referred to as “dynamic homogenization,” has attracted significant attention in recent years [246–249].

#### Applications.

As we discuss above, acoustic and elastic metamaterials draw their power from their low frequency, subwavelength properties. A nonexhaustive list of applications includes vibration and acoustic shielding, acoustic absorption, subwavelength focusing for fine-scale imaging applications, and acoustic/elastic cloaking. Below we briefly review some of the applications that have been studied. Reference [250] also provides an overview of state-of-the-art acoustic metamaterial applications.

An obvious application of a periodic material with a subwavelength bandgap is acoustic and vibration isolation, especially if this can be done over a relatively broad band of frequencies [251–253]. One may anticipate this concept finding application to sound proofing in buildings, vibration minimization in vehicle structures and home appliances, shock mitigation in armored gear, etc. In addition to isolation, where reflection is desired, researchers also considered sound absorption as an objective that also requires the minimization of reflection, e.g., see Refs. [212,254,255]. Such functionalities, would be even more effective if they can be realized in a tunable manner, which is in line with the goal of active control of phononic crystals. This problem has been pursued recently by means of fluid cavities separated by piezoelectric boundaries [256,257] and plates with a periodic array of piezoelectric transducers [258].

Another major category of applications of acoustic/elastic metamaterials is subwavelength focusing. The key principle here is utilization of operating points in the subwavelength portion of the medium's frequency spectrum whereby the isofrequency contour lines have opposing convexities across a straight path in the Brillouin zone (BZ) [208,259–265]. This concept provides a powerful alternative to conventional imaging and nondestructive inspection technologies because it essentially breaks the conventional diffraction limit.

Perhaps, the most intriguing of applications of acoustic/elastic metamaterials is cloaking, whereby an acoustic (or an elastic) wave is steered around an object rather than be scattered by it. A common approach for this problem is based on transformation acoustics (or elasticity) theory, whereby a target material distribution around the object of interest is mathematically synthesized from the point of view of heterogeneity and anisotropy [238,266–274]. Phononic crystals and acoustic/elastic metamaterials are viewed as promising candidates for practically realizing the unusually extreme material properties that emerge from this transformation process [239,271,275–284]. This is an open area of research that promises to continue to attract the attention and imagination not only of scientists, but also of the general public as a whole. For further reference, the reader may consult a review article [285] and a journal special issue [286] dedicated to the subject.

## Basic Concepts

### 1D Monatomic Lattice.

Consider the discretized elastic rod of Fig. 1 in the form of a sequence of masses connected by linear springs. A length scale d is introduced to define the distance between the masses, such that the spring constant k and mass m are, respectively, given by $k=EA/d,m=ρAd$, where E, ρ denote the elastic modulus and the density of the rod material, respectively, while A is the cross section of the rod. The position of the nth mass within the chain is expressed as xn = nd. The dynamic behavior of the resulting spring-mass chain, in the absence of an external load, is investigated by considering the equation of motion for the nth mass
$mu··n+2kun-k(un-1+un+1)=0$
(1)
where the following notation is adopted u(xn, t) = un(t). A solution is sought in the form of a harmonic plane wave
$u(xn,t)=un(t)=u∧n(ω)e-iωt=u˜(κ(ω))ei(κxn-ωt)$
(2)
where ω is the temporal frequency of harmonic motion and κ is the wavenumber, here explicitly denoted as a function of frequency. In Eq. (2), the spatial part of the solution un(ω) is expressed as
$u∧n(ω)=u˜(κ(ω))eiκxn=u˜(μ(ω))eiμn$
(3)
where $u˜$ is a complex quantity that defines the amplitude of wave motion, while the exponential term describes the phase changes as the wave propagates from one mass to the next. The amplitude $u˜$ is expressed as a function of wavenumber, and therefore of frequency, in anticipation of systems with more than a single degree of freedom in each repeated cell. In Eq. (3), the nondimensional wavenumber, also known as the “propagation constant,” μ = κd is introduced by considering the interatomic distance as the characteristic length scale of the problem, and the definition of the spatial period of the system. Equation (3) is an expression of the Floquet–Bloch theorem, which governs the propagation of plane waves in a periodic medium. A detailed discussion on the theorem can be found in Refs. [24,23]. Bloch's theorem [287] implies that the change in wave amplitude occurring from cell to cell does not depend on the cell location within the lattice, insofar as the unit cell is merely a component of an infinite lattice. Accordingly, the wave propagation characteristics of the periodic assembly can be fully identified through the analysis of the reference unit cell, provided the system does not present any dissipation, as in the case of the lattice presented here. In addition, a real wave propagation constant μ corresponds to a wave propagating without attenuation, while an imaginary component to μ defines the spatial decay in amplitude as the wave progresses through the lattice. Substituting Eq. (3) into Eq. (1) gives
$[-ω2m+2k-k(e-iμ+eiμ)]u˜(μ)ei(μn-ωt)=0$
(4)
Nontrivial ($u˜(μ)≠0$) plane-wave solutions are obtained for arbitrary locations (as defined by the mass number n) and time t by letting
$-ω2m+2k-k(e-iμ+eiμ)=-ω2m+2k(1-cosμ)=0$
(5)
Equation (5) defines the relation between frequency and wavenumber of a propagating plane wave, and therefore represents the dispersion relation for the considered 1D lattice. A nondimensional form of the dispersion relation is introduced by defining the nondimensional frequency Ω = ω/ω0, where $ω0=k/m$, which gives
$Ω2=2(1-cosμ)$
(6)

The resulting relation Ω = Ω(μ), obtained by imposing a range of wavenumbers and solving for frequency, is plotted in Fig. 2, where some of its important properties can be highlighted as discussed in the following.

Fig. 2
Fig. 2
Close modal

Periodicity. The dispersion relation is symmetric about μ = 0 and periodic with period 2π. The behavior of the relation in the fundamental period $μ∈[-π,+π]$ therefore fully describes the characteristics of harmonic plane waves in terms of the relation between their wavelength and the frequency of wave motion. This space/wavenumber duality is a feature common to all periodic media. The range of wavenumbers defining the fundamental period of the dispersion relation, i.e., $μ∈[-π,+π]$, is known as “first Brillouin zone” (FBZ) [24], which is generally defined by the reciprocal lattice, the dual of the direct lattice describing the geometry of the periodic domain. The symmetry of the dispersion relation provides the opportunity to fully characterize the dispersion properties of the lattice by only considering half of the FBZ, i.e., $μ∈[0,+π]$, which defines the “irreducible Brillouin zone” (IBZ). The periodicity also indicates that for a given value of Ω, the dispersion relation is verified for μ + 2, where m is any integer, i.e., for wavenumbers $k+2mπ/d$. The FBZ, corresponding to m = 0, is bounded by a maximum wavenumber k = π/d. This corresponds to a minimum wavelength equal to λ = 2d, which defines a wave where two neighboring masses move out-of-phase. Conceptually, all waves of wavelength λ = 2d/(2m + 1) satisfy this condition, which supports the wavelength indetermination discussed above. However, for the particular case of the discrete 1D lattice, the cases for $|m|>0$ describe wave motion with several intermass oscillations, which are not justified by the lumped parameter model and the associated degrees of freedom (DOF). This issue is, however, of great importance and introduces complexity in the analysis of periodic structures defined by unit cells with multiple degrees of freedom or of continuous periodic structures.

The periodicity of the dispersion relation can be also analyzed in light of the discretization of the rod, which essentially consists in the spatial sampling of its longitudinal displacement and physical properties at locations xn = nd. Accordingly, the sampling spatial frequency is κs = 2π/d. According to Nyquist–Shannon sampling theorem [288], this corresponds to the wavenumber range of $κ∈[-κs/2,+κs/2]$ or $κ∈[-π/d,+π/d]$, which also defines the FBZ. The Nyquist–Shannon sampling theorem states that the frequency spectrum of a sampled signal is in fact periodic with its period defined by the sampling frequency. In the case considered, this can be shown by rewriting Eq. (4) as
$[-ω2m+2k(1-cosμ)]u˜(μ)=0$
(7)
The dispersion relation for $μ'=μ+2mπ$ by virtue of its periodicity is
$-ω2m+2k(1-cosμ')=-ω2m+2k(1-cosμ)$
(8)

which leads to the following relation between wave amplitudes $u˜(μ)=u˜(μ')$, i.e., $u˜(μ)=u˜(μ+2mπ)$. Equations (2) and (3) express the plane-wave solution for harmonic propagation. In this regard, the wave amplitude is the amplitude of motion for frequency ω and wavenumber κ = μ/d, which defines the spatial spectral harmonic content of the wave at the specified frequency and wavenumber. Therefore, Eq. (8) expresses the periodicity of the spatial spectrum of the propagating wave, which is fully consistent with the spatial sampling of the continuous rod as postulated by Shannon–Nyquist theorem.

Dispersion. The nonlinearity of the dispersion relation indicates that the periodic lattice is dispersive in contrast to the continuous rod, which is nondispersive. This is another effect of the discretization process. In the long wavelength limit, the dispersion relation behaves linearly and matches that of the continuous rod, where Ω ≈ μ (Fig. 3). The phase velocity in the lattice is far from being constant and is instead given by

Fig. 3
Fig. 3
Close modal
$cph=ωk=ω0|sinμ/2||μ/2|d$
(9)

Comparison of the dispersion relations for the continuous rod and the lattice, showing the matching long wavelength behavior and of the variation of phase velocity in terms of frequency, is shown in Fig. 3. The obtained expression for the phase velocity demonstrates how the attempt to use a system of the kind of Fig. 1 to estimate the velocity of sound, leads to a complex result, while the application of partial differential calculus to the continuous system provides a much simpler value for the wave speed.

Propagation and attenuation bands. The dispersion relation solved in terms of frequency for assigned wavenumber provides solutions in a limited range, i.e., for Ω ∈ [0, 2]. For Ω > 2, there are no corresponding real-valued wavenumbers, and harmonic plane waves do not propagate. The range of Ω > 2 therefore corresponds to an attenuation range, also denoted in the literature as a “stop band” or, when this phenomenon occurs over a limited frequency range, a “bandgap.” In contrast, for Ω ∈ [0, 2], plane-wave propagation occurs without attenuation. The amount of attenuation encountered within the stop band can be quantified by solving the dispersion relation at a given frequency in terms of the wavenumber. This allows the estimation of both the real and imaginary parts of the wavenumber as functions of frequency. Figure 4 shows the imaginary part of the wavenumber, which is equal to zero in the propagation range and different than zero in the stop band. In this case, the imaginary part quantifying wave attenuation increases with frequency, which indicates that the rate of spatial decay for a wave propagating through the lattice also increases with frequency.

Fig. 4
Fig. 4
Close modal

The plots of Fig. 5 show the frequency-domain variation of the steady-state response of the last of N = 10 masses of a finite lattice excited by a harmonic force, f0 applied to the first mass. A series of resonant peaks (10) within the pass band, and strong response attenuation occurring within the stop band, are noteworthy features of the frequency response plot.

Fig. 5
Fig. 5
Close modal

The spatial attenuation is further demonstrated by plotting the dynamic deformed shapes at frequencies just below and right above the stop band boundary at Ω = 2, for a lattice of N = 100 masses. Within the pass band all masses undergo motion, which is an indication that waves propagate through the chain, while within the stop band, the amplitude of motion decays through the lattice. The magnitude of such decay, as discussed above, increases with frequency as clearly demonstrated by the deformed dynamic shapes of Fig. 6.

Fig. 6
Fig. 6
Close modal
Natural frequencies of finite systems. The dynamics of a finite lattice resulting from the “truncation” of an infinite lattice of the kind of Fig. 1 leads to a finite structure for which one can define natural frequencies and mode shapes. An estimation of the natural frequencies for a finite periodic structure can be obtained by considering the effects of truncation, or “windowing,” of the spatial extension of the lattice. This leads to the discretization of its spectral content [288], which along with the sampling effect described above, provides a good methodology for locating the natural frequencies of the chain. Due to truncation, the spectral content for the finite lattice undergoes the following transformation $u˜(μ)→u˜(nΔμ)$ or $u˜(κ)→u˜(nΔκ)$, where Δκ defines the discretization of the wavenumber domain due to the transition from the infinite to the finite lattice, and thus defines the wavenumber resolution. Such resolution is equal to
$Δκ=2πNd$
(10)
where Nd is the length of the lattice. Therefore, the natural frequencies can be read on the dispersion relation in the FBZ at values of frequencies that correspond to nΔκ, where n = 0 is included in a free–free configuration. Accordingly, the nth natural frequency is given by
$Ωn=[2(1-cos(ndΔκ))]1/2=[2(1-cos(2πnN))]1/2$
(11)
This is illustrated in Fig. 7
Fig. 7

Location on the dispersion relation of the natural frequencies of a finite, free–free lattice with N = 4 (“+”), and N = 10 (“o”) masses

Fig. 7

Location on the dispersion relation of the natural frequencies of a finite, free–free lattice with N = 4 (“+”), and N = 10 (“o”) masses

Close modal
, where the natural frequencies of a free–free lattice with N = 4 and N = 10 masses are plotted on the dispersion curve restricted to the IBZ. Further details and comprehensive analyses of the modal content of finite period structures can be found, for example, in Refs. [63,54].

### 1D Diatomic Lattice.

A simple extension to the spring-mass lattice problem is the diatomic lattice of Fig. 8, which features two different masses m1 and m2. In addition to being an interesting academic exercise, this case has physical relevance as it describes wave propagation in the [111] direction for NaCl and other crystal lattices [23].

Fig. 8
Fig. 8
Close modal
The analysis of wave motion relies on the expressions of the equations of motion for two neighboring masses (2n and 2n + 1), in the absence of external forcing and under the assumption of harmonic motion
$(-ω2m2+2k)u2n-k(u2n-1+u2n+1)=0(-ω2m1+2k)u2n+1-k(u2n+u2n+2)=0$
(12)
Equation (12) can be conveniently expressed in the following matrix form:
$(Kn-ω2M)un+Kn-1un-1+Kn+1un+1=0$
(13)

where $un=[u2n,u2n+1]T$ contains the degrees of freedom of the nth unit cell of the lattice. Under the assumption of lattice periodicity, the notation can be simplified by letting Kn = K0, where Km, with m = −1, 0, +1, are matrices defining the spring interactions within a representative unit cell and with its neighbors. Because of periodicity, these matrices are identical for all cells and are thus independent of the cell location within the lattice. Of note is the fact that the description of the system in terms of its mass and stiffness matrices suggests the extension of the techniques to more complex systems, with multiple degrees of freedom resulting, for example, from a FE discretization.

Two different approaches can be followed in order to estimate the dispersion relation of the media. The first technique generalizes the approach presented for the spring-mass system presented above, whereby the propagation constant μ is imposed to solve the resulting expression in terms of frequency. This, under the fundamental assumption of absence of dissipation, in general, leads to a simple, linear eigenvalue problem which has a convenient solution. The result of this procedure, however, does not provide any information regarding attenuation, and only predicts the real dispersion branches. This technique is often denoted as the “inverse method.” Extensions of the techniques to the analysis of systems with dissipation are an active area of research, as described in upcoming sections (see Sec. 4.1). A second technique, known as the “direct method,” leads to the dispersion relations by imposing frequency and solving for the wavenumber. The most popular implementation of the direct method formulates the “transfer matrix” of a unit cell, and evaluates its eigenvalues. The direct method is easy to implement for general 1D periodic system and can also be applied to FE discretizations of complex systems. Its main drawback consists of its difficult extension to the analysis of 2D systems, which is still an open research problem. The two techniques will be briefly described in the following.

The inverse method. The derivation of the dispersion properties for the system of Eq. (13) consists of imposing a plane-wave solution of the kind
$un(ω)=u˜(μ)einμ$
(14)
for an assigned value of the propagation constant μ. Substituting in Eq. (13) gives
$[K(μ)-ω2M]u˜(μ)einμ=0$
(15)
where $K(μ)=∑m=-1,0,+1(eimμKm)$. Equation (15) for any value of n represents a linear eigenvalue problem in terms of ω2 that can be solved for assigned values of μ in the FBZ, i.e., for $μ∈[-π,+π]$, or given symmetry, in the IBZ. The eigenvalues correspond to the various dispersion branches, while the eigenvectors $u˜(μ)$, or wave modes, describe the relative motion of the unit-cell degrees of freedom, or in general the polarization of the various wave modes. The size of the eigenvalue problem and the number of modes are determined by the number of degrees of freedom of the unit cell, which in this case is equal to 2. The two dispersion branches for the considered lattice are shown in Fig. 9,
Fig. 9

Dispersion curves of the 1D diatomic spring-mass lattice (m2 = 2m1, m1 = 1, k = 1). The lower branch represents acoustic modes (solid line) and the upper branch represents optical modes (dashed line).

Fig. 9

Dispersion curves of the 1D diatomic spring-mass lattice (m2 = 2m1, m1 = 1, k = 1). The lower branch represents acoustic modes (solid line) and the upper branch represents optical modes (dashed line).

Close modal
, for the case of m2 = 2m1, m1 = 1, k = 1. The lower branch, which in the long wavelength coincides with the dispersion branch of an acoustic wave in a homogeneous medium, is known as the “acoustic mode,” while the top branch is called the “optical mode.” Of note is the fact that this lattice is again characterized by a stop band above the optical mode and by a bandgap that separates the acoustical and optical mode, whose width is defined by the mass ratio m2/m1, and becomes equal to zero in the limiting case of m1 = m2, which reduces to the lattice of Fig. 1.

The direct method and transfer matrix. The implementation of the direct method generally relies on the formulation of the transfer matrix of the lattice, which relates the generalized displacements and forces at the left and right cell boundaries (see Fig. 10). The description of the behavior of the unit cell in Fig. 10 in terms of the dynamic stiffness matrix $D(ω)=K-ω2M$ is a convenient starting point to relate boundary forces to boundary displacements. In expanded form, the dynamic stiffness relation can be written as

Fig. 10
Fig. 10
Close modal
${fLfR}n=[DLLDLRDRLDRR]{uLuR}n$
(16)
where $DRL=DLRT$. This relation can be rearranged to relate the generalized state vectors, containing displacement and forces, at the cell left (L) and right (R) boundaries. This gives
${uRfR}n=[-DLR-1DLLDLR-1DRRDLR-1DLL+DRLDRRDLR-1]{uLfL}n$
(17)
Next, the state vector at the right boundary of cell n can be related to the state vector at the left boundary of cell n + 1 by imposing displacement compatibility and force equilibrium at the interface, i.e., $uRn=uLn+1,fRn=-fLn+1$. The state vectors at the same location (say the left boundary) in two consecutive cells are therefore related as follows:
$yn+1=Tyn$
(18)
where y = [u, f]T is the state vector, and
$T=[-DLR-1DLLDLR-1DRRDLR-1DLL-DRL-DRRDLR-1]$
(19)

is the transfer matrix.

Wave propagation through the lattice can be investigated by applying Floquet's theorem, which seeks solutions where the state vectors of two consecutive cells are related by
$yn+1=λyn$
(20)
where λ is a Floquet multiplier, which is related to the propagation constant by letting λ = e. Comparing Eqs. (19) and (20) leads to the following eigenvalue problem:
$Tyn=λyn$
(21)
whose solution provides the values of the Floquet multipliers at the considered frequency. From the Floquet multipliers, the propagation constants are easily determined so that their real and imaginary parts are evaluated as functions of frequency. Results for the considered diatomic spring-mass lattice are presented in Fig. 11
Fig. 11

Dispersion curves for the 1D diatomic spring-mass lattice (m2 = 2m1, m1 = 1, k = 1) showing real (solid line) and imaginary (dashed line) parts of the propagation constant

Fig. 11

Dispersion curves for the 1D diatomic spring-mass lattice (m2 = 2m1, m1 = 1, k = 1) showing real (solid line) and imaginary (dashed line) parts of the propagation constant

Close modal
, where acoustic and optical branches are plotted along with the corresponding attenuation curves. The approach presented herein for the derivation of the transfer matrix is general and can be applied to systems discretized through a FE procedure, as well as to continuous systems whose transfer matrix can be formulated on the basis of analytical equations of motion [289,290]. The solution of Eq. (21) can leverage several of the interesting properties of the transfer matrix, which are reviewed, for example, in Refs. [291–295].

### 1D Lattice With Internal Resonator.

A variant to the case of the 1D monatomic lattice of Fig. 1 consists of an addition of a second-order oscillator of spring constant kR and mass mR, connected to each base mass (Fig. 12). This leads to a configuration characterized by internal resonances centered at the resonant frequency of the oscillators, i.e., $ωR=kR/mR$. This produces a bandgap in the vicinity of ωR. In contrast to the bandgaps discussed for monatomic and standard diatomic lattices, the internal resonance gap is not related to wavelength, and therefore in principle can be produced at arbitrarily low-frequencies/long wavelengths. The bandgaps observed in monatomic and standard diatomic lattices instead only occur at wavelengths which are of the order of the unit-cell size, and therefore are strictly constrained by the unit-cell dimensions. These are the result of interference generated by the interaction of incident and scattered waves at the unit-cell boundaries, which is the phenomenon often referred to as Bragg scattering. The wavelength-independent nature of internal resonance bandgaps has important implications, which have motivated several studies for practical implementation of the concept [22], along with numerous studies devoted to the interpretation of the physical meaning of the locally negative slope of the dispersion relations caused by the internal resonance, which is an indication of negative group velocity [224]. Interpretations in terms of negative mass can be found, for example, in Refs. [240,243,296]. The reader is referred to Sec. 1.4.3 for a detailed overview of the literature on locally resonant metamaterials, and to other parts of this review where the notion of local resonance has been considered.

Fig. 12
Fig. 12
Close modal
The analysis of the dispersion properties of the 1D chain in Fig. 12 is conducted by considering the motion of unit cell n, which is governed by the following equations of motion:
$(-ω2m+2k)un-k(un-1+un+1)-kR(uRn-un)=0(-ω2mR+kR)uRn-kRun=0$
(22)
where mR, kR, respectively, are the mass and spring constant of the oscillator, while uRn is the associated degree of freedom. Equation (22) can be combined by condensing the resonator degree of freedom, which gives
$(-ω2m+2k-kR2kR-ω2mR+kR)un-k(un-1+un+1)=0$
(23)
Application of Bloch's theorem and a few manipulations lead to the following dispersion relation:
$2(1-cos(μ))-Ω2(1+m¯R1-Ω2/ΩR2)=0$
(24)
where Ω = ω/ω0 (as previously defined) and ΩR = ωR/ω0, while $m¯R=mR/m$ defines the ratio of the mass of the resonator to that associated with the primary chain. Solution of the dispersion relation above for harmonic wave motion, i.e., for assigned values of frequency, provides the real and imaginary parts of the wavenumber μ. The presence of the internal resonators causes the dispersion curve to split into two separate branches, between which a bandgap is observed. This bandgap is centered at the tuning frequency ΩR, over a range which is also highlighted by the nonzero attenuation parameter (imaginary part of the propagation constant). The presence of dissipation in the system will affect the dispersion diagram and will cause the two branches to merge into a single curve featuring a characteristic “S” shape in a frequency range centered at the internal resonance. The corresponding backward slope of the dispersion branch defines a range where the group velocity is negative, and is one distinctive feature of internally resonating metamaterials, of which the lattice considered herein is an elementary example. Figure 13,
Fig. 13

Dispersion curves for a 1D spring-mass lattice with internal resonators. Comparison between undamped curves obtained for $m¯r=0.125,Ω0=1,ΩR=1$ (black dashed line), and “S” shaped branches corresponding to a viscously damped configuration (red solid line).

Fig. 13

Dispersion curves for a 1D spring-mass lattice with internal resonators. Comparison between undamped curves obtained for $m¯r=0.125,Ω0=1,ΩR=1$ (black dashed line), and “S” shaped branches corresponding to a viscously damped configuration (red solid line).

Close modal
compares the dispersion curves obtained for the undamped system as expressed by Eq. (24) (black dashed line) with that of a damped system with linear viscous dampers connecting the masses (red solid line) (see upcoming Fig. 24). The figure illustrates the split of the dispersion branches in the undamped case, and the “S” shape of the curve associated with the damped case. Additional details on similar damped configurations and on the effects of damping on free wave motion are provided in Sec. 4.1.3.

The frequency of resonance ΩR can be selected, or “tuned,” to obtain attenuation over a desired band, without any constraint imposed by wavelength. Of course, in practice, low-frequency attenuation requires large mass values and/or low stiffness, which poses a challenge for the implementation of this concept. However, the application of the internal resonator concept for systems more complex than the one illustrated in this section has been broadly pursued by the research community, and several experimental demonstrations are also available. Examples of configurations and experimental results are presented in Sec. 5 of this paper.

### 2D Monatomic Lattice.

The spring-mass lattice can be extended to obtain a 2D periodic structure that develops over a plane (Fig. 14). The unit cell consists of a mass undergoing motion perpendicular to the plane of the structure and connected to its neighbors through two springs that apply a force proportional to the relative out-of-plane motion of two neighboring masses. This example, which was originally introduced in Ref. [76], can be described through a single degree-of-freedom unit cell, which is sufficient to uncover additional interesting properties of periodic systems. According to the configuration of Fig. 14, wave motion is investigated by expressing the equation of a generic unit cell at location m, n

Fig. 14
Fig. 14
Close modal
$mu··m,n+2(kx+ky)um,n-kx(um,n-1+um,n+1) -ky(um-1,n+um+1,n)=0$
(25)
Similarly to the 1D case, a plane-wave solution is sought in the form
$um,n(t)=u∧m,n(ω)e-iωt=u˜(κ(ω))ei(κ·xm,n-ωt)$
(26)
where $κ=κxi+κyj$ is the wavevector, while $xm,n=maxi+nayj$ defines the position of the unit cell in the lattice. Two length scales ax, ay are introduced to define the interatomic distances along the two in-plane directions. The spatial part of the solution um,n(ω) can be rewritten as
$u∧m,n(ω)=u˜(μ(ω))ei(mμx+nμy)$
(27)
where $μ=μxi+μyj$ is the propagation vector, whose components are $μx=κxax,μy=κyay$. Substituting Eq. (27) into Eq. (26) gives
$[-ω2m+2kx(1-cosμx)+2ky(1-cosμy)]u˜(μ)=0$
(28)
upon following steps similar to those used for the 1D spring lattice. The existence of nontrivial harmonic plane-wave solutions requires that
$-ω2m+2kx(1-cosμx)+2ky(1-cosμy)=0$
(29)
which defines the dispersion relation for the considered 2D lattice. Assuming for simplicity that kx = ky = k allows rewriting Eq. (29) in the following nondimensional form:
$Ω2=2[2-cosμx-cosμy]$
(30)
where Ω = ω/ω0 is the nondimensional frequency. The nondimensional dispersion relation describes a surface in terms of the components of the propagation vector, i.e., $Ω=Ω(μ)=Ω(μx,μy)$. Figure 15,
Fig. 15

Dispersion surface of the 2D spring-mass lattice with kx = ky and ax = ay

Fig. 15

Dispersion surface of the 2D spring-mass lattice with kx = ky and ax = ay

Close modal
shows a color map of the dispersion surface, where the third dimension is frequency Ω. The plot clearly illustrates the periodicity of the surface in the wavenumber domain, and highlights the FBZ and the IBZ. The fundamental period defines the FBZ, which is defined by $μx,μy∈[-π,+π]$. In addition, the symmetry of the surface can be exploited to identify an IBZ, whose edges are identified by the letters Γ, X, M in Fig. 15. Compact representations of the dispersion properties of multidimensional periodic systems are obtained by plotting frequency as the wave vector's components vary along the boundary of the IBZ, which leads to a widely accepted and effective visualization of bandgaps and overall dispersion properties.
In addition to bandgaps, the dispersion surfaces provide information on the propagation of energy within the lattice, which is defined by the group velocity
$cg=∇κω$
(31)
For the 2D lattice considered herein, one can conveniently show that the directions of energy flow within the structure at a given frequency are identified by the normals to the associated isofrequency contour line of the dispersion surface. The components of the group velocity $cg=cgxi+cgyj$ for the considered 2D lattice, assuming ax = ay = a and kx = ky = k, are given by
$cgx=a∂ω∂μx=ω0a2sinμx[2-cosμx-cosμy]1/2cgy=a∂ω∂μy=ω0a2sinμy[2-cosμx-cosμy]1/2$
Hence, the direction of the group velocity vector is given by
$tanφ=sinμysinμx$
(32)
The considered isofrequency line can be expressed in the following functional form μy = f(μx), which describes a curve in the μx, μy plane. Its local slope at a given frequency is given by
$∂μy∂μx=-sinμxsinμy=tanψ$
(33)
Thus, $tanψ=-(1/tanφ)$ or $ψ-φ=(π/2)$, which indicates that the local slope of the dispersion surface, as defined by the direction of the local normal to the isofrequency line, defines the direction of wave propagation at the considered wavenumber pair μx, μy, and frequency ω. The identification of the direction of energy propagation through the analysis of the slopes of the dispersion surface is illustrated in Fig. 16,
Fig. 16

Dispersion surface contours and arrows showing the normal directions to isofrequency contours at Ω = 1 and $Ω=2$

Fig. 16

Dispersion surface contours and arrows showing the normal directions to isofrequency contours at Ω = 1 and $Ω=2$

Close modal
. The dispersion surface is represented as a contour plot, and the directions of wave motion at two selected frequencies (Ω = 1 and $Ω=2$) are highlighted by arrows normal to the corresponding isofrequency contours. The Ω = 1 isofrequency contour appears as almost circular, which is indicative of the fact that energy at this frequency is expected to propagate in all directions almost uniformly. In fact, a perfectly circular isofrequency contour line would correspond to a perfectly isotropic behavior at that frequency. The case for $Ω=2$ is remarkably different as the isofrequency contours are composed of straight lines, and all normal arrows are directed solely at ±45 deg. This indicates that at this frequency, the lattice is characterized by strong anisotropy and directivity of wave motion which is confined to the ±45 deg directions. This is confirmed by the plots of the group velocity at the two frequencies, which is almost circular at Ω = 1, while it features a star shape with strongly focused lobes at $Ω=2$ at ±45 deg (Fig. 17
Fig. 17

Variation of group velocity at Ω = 1 (black dotted line) and $Ω=2$ (red solid line)

Fig. 17

Variation of group velocity at Ω = 1 (black dotted line) and $Ω=2$ (red solid line)

Close modal
).

The strong directivity is further illustrated by computing the harmonic response of a finite lattice at the chosen frequencies. A finite 81 × 81 lattice in free–free configuration is excited by a concentrated load applied at the center of the lattice. The steady-state harmonic response is computed by inverting the dynamic stiffness matrix of the assembly, and is plotted in the form of the color map of the response magnitude. Proportional damping is added in order to reduce the effects of boundary reflections and better visualize the response directivity. The response plots for the two considered frequency values are shown in Fig. 18, which demonstrate the difference in dynamic behavior, and, in particular, the strongly directional response of the lattice at $Ω=2$.

Fig. 18
Fig. 18
Close modal

The occurrence of wave “beaming” due to limited directionality associated with specific frequencies of wave propagation is illustrated in this simple case, and has been broadly investigated in the literature. Directional wave propagation in crystal structures of various kinds of anisotropy has been studied extensively and documented both numerically and experimentally [76]. Strong focusing of waves in specific directions is associated with the occurrence of “caustics” [84,297], which correspond to closed circles in the group velocity plots as shown in Refs. [80,298,299]. Directional wave propagation has been investigated in periodic structural components, such as beam grillages, hexagonal, re-entrant and chiral lattices, and periodic plates among others. In addition, the analysis of dispersion surfaces and group velocity directions has provided the foundation for the design and implementation of phononic crystals with negative acoustic refraction and focusing capabilities [156,160,163,260,261,300, 301,302].

## Computational Analysis Techniques

While modeling macroscale periodic media as lumped parameter spring-mass systems provides a simplification and hence insight into the dynamical characteristics, it is often more efficient and accurate to employ continuum modeling. The governing equation of motion for a linear phononic material modeled as a continuum is
$∇·C:∇Su=ρu··$
(34)
where u is the displacement, C is the elasticity tensor, ρ is the density, and ∇S is the symmetric gradient operator, i.e.,
$∇Su=12(∇u+(∇u)T)$
(35)
Applying Eq. (34) to a unit cell representing a phononic crystal (or an acoustic elastic metamaterial) and assuming a Bloch wave solution of the form
$u(x,κ;t)=u˜(x,κ)ei(κ·x-ωt)$
(36)

allows us to obtain the frequency band structure along any desired path in, or at the boundaries of, the BZ. In Eq. (36), $u˜$ is the unit-cell Bloch displacement function that is periodic over the unit cell, x = {x, y, z} is the position vector, and $κ={κx,κy,κz}$ is the wavevector.

### Overview of Band Structure Computational Techniques.

This section reviews a variety of computational techniques for analyzing response and designing continuous phononic systems. Included are descriptions of the plane-wave expansion methods (PWE), finite difference (FD) methods, the FE methods, and the multiple scattering (MS) method. The literature reviewed is not intended to list exhaustively those papers using each method, but instead to document early literature employing the method and extending it, and to provide the proper technical review necessary for a reader to begin using each method.

#### Plane-Wave Expansion Method.

The general approach used in the plane-wave expansion method for continuous systems is to expand the sought-after solution field and the material properties in a Fourier series, and to then invoke orthogonality of the basis functions to solve individually each of the introduced solution coefficients. A very early description of this method can be found in the beginning of Sec. 2 of Elachi's review article [85], although the method is not given a name therein. Phononic crystal research in which this approach has been taken date back to the early 1990s [17,18] with work on acoustic and elastic composites.

The basic PWE ideas can be illustrated for a one-dimensional, continuous system. Consider such a system governed by the periodic wave equation
$ρ∂2u∂t2=∂∂x(ρc2∂u∂x)$
(37)
where u(x, t) denotes the displacement field, ρ(x) the density, and c(x) the speed of sound through the media. The material properties ρ(x) and c(x) are assumed to vary periodically with periodicity a, i.e., the lattice constant. A lattice vector for this system can be written as $R=na=nax∧$ where $x∧$ denotes a unit vector and n any integer. Similarly, a reciprocal lattice vector can be written as $Gm=mb=m/ax∧$ where m denotes any integer. The PWE method introduces expansions
$u(x,t)=ei(κ·r-ωt)∑Gluκ(Gl)eiGl·r$
(38)
$ρ(x)=∑Gmρ(Gm)eiGm·r$
(39)
$ρ(x)c2(x)=∑Gmτ(Gm)eiGm·r$
(40)

where $r=xx∧$ denotes position, $κ=kx∧$ denotes the wavevector, and $∑Gl,∑Gm$ denote sums over all reciprocal lattice vectors (e.g., over all l, m). The transformed density and stiffness are denoted by $ρ(Gm)$ and $τ(Gm)$, respectively. The transformed displacement is given by $uκ(Gl)$ for any wavevector $κ$. Note that Eq. (38) embodies Bloch's theorem.

Substituting the expansions into Eq. (37) and forming the complex inner product with $eiGn·r$ (its complex conjugate appearing in the integral) yields
$∑Gl[-ω2ρ(Gn-Gl)+(κ+Gl)·(κ+Gn)τ(Gn-Gl)]uκ(Gl)=0$
(41)

Used in arriving at Eq. (41) is the fact that for each $Gn$, the only nonzero terms in the inner product are those satisfying $Gl+Gm-Gn=0$, leading to the condensing of the summation over $Gm$. By truncating the expansions for $Gl$, and evaluating Eq. (41) for a similarly truncated set of $Gn$, a well-posed eigenvalue problem results for eigenfrequencies $ω(κ)$ and eigenvectors $[uκ(G-N)uκ(G-N+1)…uκ(GN)]T$, where 2N + 1 denotes the number of retained terms. The band structure is then obtained by sampling $κ$ in the first Brillouin zone. As in any approximation technique, a convergence criterion should be used to determine an appropriate value for N. Convergence is discussed in Ref. [303], where it is pointed out that it may be appropriate to expand the inverse of $ρ(x)c2(x)$ to achieve faster convergence.

#### Finite Difference Methods.

The finite difference method is a classical method that has seen wide application in almost every discipline where partial differential equations (PDEs) govern behavior. It has one critical advantage over PWE methods in that it can easily analyze periodic media composed of multiple states (e.g., solid and fluid) [304]. This includes a large class of materials with voids. The earliest use in a phononic system known to the authors is that by Yang and Lee [305]. They found the band structure of a two-layered one-dimensional composite (what is now referred to as a bimaterial rod) assuming harmonic time dependence and using a spatial finite difference approximation of the resulting governing equation. The authors also enforce a Floquet (or Bloch) condition on the meshed domain, effectively limiting the computation to a single unit cell. Solution of a discrete eigenvalue problem at each frequency yields the dispersion curves.

A more frequently encountered finite difference approach is one in which a PDE is discretized in both space and time to compute the time-domain response. Frequency content is then acquired from Fourier transforms (FTs) of the time-history data. This subapproach is appropriately termed the finite difference time-domain (FDTD) approach. Note that if the Floquet (or Bloch) condition is not imposed in arriving at the PDE, the computational domain must include many unit cells, an excitation source, and a choice of boundary conditions. For two-dimensional phononic crystals, such an approach was first presented by Garcia-Pablos et al. [304]. By varying the source's excitation frequency, attenuation versus frequency can be characterized and used to identify bandgaps. In contrast, Tanaka et al. [306] presented an FDTD approach in which the PDE is formulated with the Floquet (or Bloch) condition, thus requiring simulation only over a unit cell. By varying the wavevectors through the Brillouin zone, Fourier transforms of the time-history data can be used to generate band structures.

The use of the frequency-domain finite difference method for generating band structure is overviewed next for the periodic wave equation Eq. (37). Assuming time-dependence eiωt and introducing centered difference approximations $(∂2u/∂x2)≈(ul+1-2ul+ul-1)/h2$, where h denotes the discretization length and l denotes a grid point, the differenced system of equations for a single unit cell can be written as
$-ω2AU=1h2BU$
(42)
where U holds the displacements ul, $l=0…n+1$; n denotes the number of grid points unique to the unit cell; and A and B are diagonal and tridiagonal matrices, respectively, containing material properties associated with each grid point. Frequently, the stiffness assigned to a grid point is that at the forward half grid point, and the density is that at the grid point itself, where care is taken not to allow a grid point to lie directly on an interface [305]. Other choices are, however, possible. Note that U holds two grid points assigned to neighboring unit cells, namely, u0 and un+1. These are removed with the Bloch conditions
$u0=e-iμun,un+1=eiμu1$
(43)

Substitution of Eq. (43) into Eq. (42) yields the desired n by n complex eigenvalue problem. Sweeping through μ = ka and solving the eigenvalue problem at each wavenumber generates the band structure.

#### Finite Element Method.

The finite element method is also a classical method for determining solutions to PDEs. With respect to phononic systems, the first study to employ it in generating band structures appears to be that by Orris and Petyt [70]. The linear, undamped version of the method yields a set of matrix equations for a single unit cell that can be posed as
$-ω2Mq+Kq=0$
(44)

where M denotes a mass matrix, K denotes the stiffness matrix, and q stores the nodal degrees of freedom. Note that these equations are arrived at by considering a unit cell in isolation from other unit cells, and thus neighbor forces do not appear on the right-hand side; other treatments can be found in Refs. [82,307,308] where neighbor forces are included, but lead to the same final form. As can be expected, the finite difference equations Eq. (42) take this same form if rearranged with the discretization length h incorporated into the stiffness matrix; one notable difference is that in the FE treatment, the mass matrix is not necessarily diagonal. Since the one-dimensional case was discussed for finite difference modeling, the two-dimensional case will be discussed here for the finite element case.

In arriving at Eq. (44), all elements lie wholly within the unit cell and all boundary nodes lie on the perimeter of the unit cell (i.e., not inside of the unit cell). Thus, all unit-cell energy is accounted for, and no energy associated with M or K belongs to a neighboring unit cell. However, an excess number of nodes exist in q. Specifically, nodal degrees of freedom at the top and left sides of the unit cell (see Fig. 19) need to be condensed out since they belong to neighboring unit cells.

Fig. 19
Fig. 19
Close modal
For plane waves, the Bloch conditions reduce the nodal degrees of freedom
$q=Tq∧$
(45)
where T denotes a linear transformation parameterized by the wavenumbers μx and μy. For example, denoting the reduced set of nodal displacements by $q∧$, we have $q=Tq∧$ where
$q=[qiqBqTqLqRqLBqRBqLTqRT], T=[I0000I000Ieμy0000I000Ieμx0000I000Ieμx000Ieμy000Ieμx+μy], q∧=[qiqBqLqLB]$
(46)
Substitution of Eq. (46) into Eq. (44) results in an equation of motion in the form of $(-ω2M+K)Tq∧=0$. Introducing $T¯$ obtained from T by replacing μx, μy with −μx, −μy, and premultiplying the equation of motion by $T¯T$, the final equations of motion are formed
$(-ω2M∧+K∧)[qiqBqLqLB]=0$
(47)

where $M∧(μx,μy)=T¯TMT$ and $K∧(μx,μy)=T¯TKT$. Sweeping through all values of μx and μy along the irreducible Brillouin zone and solving the eigenvalue problem generate the band structure.

An alternative FE approach is based on undertaking a Bloch operator transformation of the governing differential equations, which is done by substituting Eq. (36) into Eq. (34) to obtain the strong form of the Bloch eigenvalue problem [309]
$∇·C:[∇Su˜+i2(κT⊗u˜+κ⊗u˜T)]=-ρω2u˜$
(48)
for which standard periodic boundary conditions are applied. We note that in Eq. (48), the wavevector is denoted by $κ$. Using variational principles, the weak form for Eq. (48) is
$∫Ω[∇Sw˜+i2(κT⊗w˜+κ⊗w˜T)]:C: [∇Su˜+i2(κT⊗u˜+κ⊗u˜T)]dΩ=ω2∫Ωρ¯w˜·u˜dΩ$
(49)
where $w˜$ is a Bloch-transformed weighing function. The operator on the left-hand side of Eq. (49) arising from the transformation is Hermitian. In general, the right-hand side operator may also be Hermitian, as this depends on the form of $w˜$. Equation (49) may now be discretized using the FE method and upon application of regular periodic boundary conditions the unit-cell Bloch eigenvalue problem is obtained
$(K'(κ)-ω2M')U˜=0$
(50)

where $U˜$ is the discrete Bloch vector. This equation may be solved for different values of $κ$ to generate the band structure. More details on the FE discretization of Eq. (49) are available in Ref. [309].

#### Multiple Scattering Method.

The MS method appears to be first proposed for evaluating band structure in periodic elastic composites by Economou and coworkers starting in 1999 [310,311] and by Liu and coworkers in 2000 [312]. As the name of the method implies, it is particularly useful in phononic systems composed of scatterers in a host medium (in contrast to layered phononic systems). The method can handle multiple media states (e.g., elastic scatters in a fluid, or air spaces in an elastic medium) and can accurately analyze high-contrast problems where the other methods (PWE, FD, FE) may exhibit slow convergence.

An excellent review of the multiple scattering method and its use in determining both band structure and frequency-dependent transmission coefficients can be found in Ref. [313]. Reviewing the method, even in simple cases, requires a lengthy development beyond the scope of the present review. For the purposes herein, we only present the general ideas that guide the method.

Considering the scalar case of acoustic waves scattered by fluid or solid scatterers, the method begins by decomposing the total three-dimensional pressure field p(r) into an initially incident p0(r) and scattered field psc(r), where the scattered field is composed of multiple scattered components due to scatterers at locations Rn
$p(r)=p0(r)+∑npnsc(r)$
(51)

Note that a scatterer located at Rn generates a pressure field $pnsc(r)$ that in itself arises from incident pressure waves from neighboring scatterers. If the single scatterer problem has been solved (i.e., the transmission coefficient is known for an incident spherical wave on the scatterer), a linear algebraic system can be established that determines amplitude coefficients for each scattered wave leaving a scatterer at Rn. For an acoustic medium, the propagation away from the scatterer is then the product of the amplitude coefficients with spherical Hankel and spherical harmonic functions whose arguments contain rRn, as per a standard Green's function approach. Note that the problem as described is infinite in size, so computed solutions require truncating the number of neighbors used to describe the waves scattered at Rn. Finally, if band structure is desired, Bloch's theorem can be invoked to relate amplitude coefficients for a scatterer in one unit cell to those of an analogous scatterer in a reference unit cell.

### Model Reduction for Band Structure Calculations.

As thoroughly discussed in Secs. 1 and 2, the study of wave propagation in periodic media in general utilizes Bloch's theorem, which allows for the calculation of dispersion curves (frequency band structure) and density of states. Due to crystallographic symmetry, the Bloch wave solution needs to be applied only to a single unit cell in the reciprocal lattice space covering the first BZ [24]. Further utilization of symmetry reduces the solution domain, even more, to the IBZ. As mentioned in Sec. 3.1, there are several techniques for band structure calculations for phononic crystals and acoustic/elastic metamaterials (which are also applicable to photonic crystals and electromagnetic metamaterials). Some of the methods involve expanding the periodic domain and the wave field using a truncated basis. This provides a means of classification in terms of the type of basis, e.g., the PWE method (Sec. 3.1.1) involves a Fourier basis expansion and the FE method (Sec. 3.1.3) involves a real space basis expansion. The pros and cons of the various methods are discussed in depth in the literature [314].

Regardless of the type of system and type of method used for band structure calculations, the computational effort is usually high because it involves solving a complex eigenvalue problem and doing so numerous times as the value of the wavevector $κ$ is varied. Therefore, means for speeding up band structure calculations are very beneficial. Examples of techniques that achieve this speed-up include utilization of the multigrid concept [315], development of fast iterative solvers for the Bloch eigenvalue problem [316,317], extension of homogenization methods to capture dispersion across the spectrum [43,318,44], and formulation of multiscale finite elements as illustrated recently in Ref. [319].

In this section, we describe a model reduction method that is based on modal transformation [309]. This method, which is referred to as the reduced Bloch mode expansion (RBME) method, involves carrying out an expansion employing a natural basis composed of a selected reduced set of Bloch eigenfunctions. This reduced basis is selected within the IBZ at high-symmetry points determined by the crystal structure and group theory (and possibly at additional related points). At each of these high-symmetry points, a number of Bloch eigenfunctions are selected up to the frequency range of interest. As mentioned above, it is common to initially discretize the problem at hand using some basis of choice. In this manner, RBME constitutes a secondary expansion using a set of Bloch eigenvectors and hence keeps and builds on any favorable attributes the primary expansion approach might exhibit. The proposed method is in line with the well-known concept of modal analysis, which is widely used in various fields in the physical sciences and engineering. To follow, a brief description of the RBME process is given and its application in a discrete setting (e.g., using FEs) is demonstrated for a phononic crystal problem.

The starting point for the RBME method is a discrete generalized eigenvalue problem emerging from the application of Bloch's theorem. For demonstration, Eq. (50) will be used for this purpose. The Bloch eigenvalue problem is first solved at a reduced set of selected wavevector points (i.e., reduced set of $κ$-points). This provides the eigenvectors from which a reduced Bloch modal matrix, denoted Ψ, is formed. Several schemes are available for the selection, the simplest of which is the set of eigenvectors corresponding to the first few branches at the high-symmetry points, e.g., the Γ, X, and M points for a 2D square lattice as illustrated in Fig. 20 (more details on selection schemes are given in Ref. [309]). The matrix Ψ is then used to expand the eigenvectors $U˜$, i.e.,

Fig. 20
Fig. 20
Close modal
$U˜(n×1)=Ψ(n×m)V˜(m×1)$
(52)
where $V˜$ is a vector of modal coordinates for the unit-cell Bloch mode shapes. In Eq. (52), n and m refer to the number of rows and number of columns for the matrix equation. To enable significant model reduction, the chosen $κ$-point selection scheme has to ensure that mn. Substituting Eq. (52) into Eq. (50), and premultiplying by the complex transpose of Ψ, yields a reduced eigenvalue problem of size m × m
$Ψ*K'(κ)ΨV˜-ω2Ψ*M'ΨV˜=0,KR(κ)V˜-ω2MRV˜=0$
(53)

where MR and $KR(κ)$ are reduced generalized mass and stiffness matrices. The eigenvalue problem given in Eq. (53) can then be solved for the entire region of interest within the IBZ at a significantly lower cost compared with using the full model given in Eq. (50).

To demonstrate the RBME approach, we analyze a linear elastic, isotropic, continuum model of a 2D phononic crystal under plain-strain conditions (the method is also applicable to periodic metamaterials). As an example, a square lattice is considered with a bimaterial unit cell (with ratios of Young's moduli and densities of E2/E1 = 16 and ρ2/ρ1 = 8, respectively). The topology of the material phase distribution in the unit cell is shown in the inset of Fig. 20 (the filling ratio of the inclusion is 0.1111). The unit cell is discretized into 45 × 45 uniformly sized four-node bilinear quadrilateral FEs, i.e., 2025 elements. With the application of periodic boundary conditions, the number of degrees of freedom is n = 4050. Figure 20 shows the calculated band structure (a) and density of states (b) using two-point expansion, that is, the selection is carried out at the Γ, X, M points in $κ$-space. In the calculations, eight modes were utilized at each of these selection points. As such, a total of 24 eigenvectors (m = 24) were used to form the Bloch modal matrix. The results for the full model are overlaid for comparison indicating excellent agreement, despite a reduction of model size from 4050 to 24 degrees of freedom. For models with a larger number of degrees of freedom, and a calculation with high $κ$-point sampling, two orders of magnitude or greater reduction in computational expense will be achieved (as shown in Fig. 20(c)).

## Theoretical Studies

### Damped Phononic Materials and Structures.

The presence of damping is unavoidable in materials and structures. Its consideration therefore improves the predictive ability of any model of a phononic system [320]. From a characterization perspective, damping can be identified from experiments [321]. From a design perspective, damping may be introduced intentionally to provide stronger attenuation and/or more control of the wave propagation characteristics, as will be demonstrated in this section. In both the prediction and design scenarios, the treatment of damping within the framework of dispersion analysis provides a fundamental description of the effects of dissipation on the intrinsic properties of a phononic material or structure. This is more powerful than consideration of the effects of damping on finite systems where only extrinsic (i.e., structural level) response behavior may be observed.

In this section, we consider the problem of damped phononic materials and structures. We first review relevant studies from the literature in Sec. 4.1.1, and follow in Secs. 4.1.2 through 4.1.4 with a detailed technical review on the treatment of viscous damping in the calculation of the dispersion curves. The technical review will consider (1) a 1D diatomic chain like the one depicted in Fig. 8 except with viscous damping elements included between the masses; (2) a 1D diatomic chain incorporating local resonators similar to the one depicted in Fig. 12 but with damping included; and (3) a 2D continuum model of a damped phononic crystal discretized by the finite element method.

#### Literature Review.

There are numerous studies in the literature that consider the treatment of damping in the context of periodic materials or structures. Some of these studies are concerned with finite structures [322–324] or the transmission [325] or Green function [326] response in a portion of a periodic material or structure. In this review, the interest is in the examination of materials and structures mathematically modeled as infinite systems since this provides intrinsic wave propagation properties as discussed above. Research in this area may be divided into two categories. The first seeks a wave dispersion solution for prescribed harmonic motion and the second is concerned with obtaining a dispersion relation for free wave propagation. In practice, the first problem is useful for comparison with experiments on long phononic structures, or portions of a phononic material, subjected to forced harmonic excitation. The second problem is suited for investigation of the natural dynamic characteristics of phononic materials and structures or when the loading conditions are of an impulsive nature.

A key distinction between harmonic and free wave propagation lies in whether the frequency is restricted to real values (harmonic) or allowed to be complex (free). The wavenumber (or wavevector) can be either, depending on whether the interest is in both propagating and evanescent modes (for which the wavenumber has to be complex) or in only propagating modes (for which it is sufficient to only consider a real wavenumber).

The study of harmonic wave propagation in damped periodic materials and structures in the context of a linear eigenvalue problem goes back to Mead [6] in which a monocoupled periodic chain with a complex stiffness connecting each unit was considered. Langley [327] investigated a similar model only with the damping introduced in the inertial term. Merheb et al. [328] and Zhao and Wei [329] also incorporated damping in a linear eigenvalue problem setting but considered more complex frequency-dependent damping models that carried information on the time history of the dissipation mechanism. The approach by Merheb et al. [328] was based on the FD time-domain method, whereas Zhao and Wei [329] used the PWE method. Recently, Manconi and Mace [330] used finite elements to treat frequency-dependent damping in the context of harmonic wave propagation in laminated panels. On a different track, treatment of harmonic wave propagation through the setting up of a quadratic eigenvalue problem has been considered for both 1D models [331–333] and multidimensional models [334–336]. The advantages of the quadratic eigenvalue problem route are twofold: (1) it easily enables the incorporation of frequency dependency in the material properties and (2) it readily provides the wave solution for both spatially propagating and spatially decaying modes.

Calculation of dispersion relations for free wave propagation in damped phononic materials and structures may be traced as far back as the 1975 paper by Mukharjee and Lee [337] that studied laminated media with structural damping and where the finite difference method was used to obtain dispersion relations for various levels of damping. Sprik and Wegdam [338] and Zhang et al. [339] later considered similar problems—in which damping was introduced by adding an imaginary component to the elastic modulus—using the plane-wave expansion method. Recently, Hussein [340] and Hussein and Frazier [341] derived dispersion relations for viscously damped phononic materials. The treatment of temporally nonlocal viscoelasticity was also considered [342]. In all studies concerned with free wave propagation, a complex frequency is assumed in the application of Bloch's theorem. This provides the usual frequency dispersion curves—now affected by the presence of damping—and also curves associated with the temporal attenuation for each Bloch mode. A wavenumber-dependent damping ratio diagram may be extracted from the latter to explicitly determine the loss factor for each Bloch mode [340–342].

#### Damped 1D Diatomic Lattice.

To demonstrate the treatment of damping in free wave propagation analysis, we start with the simple 1D diatomic lattice model studied in Sec. 2.2, and add linear viscous damping elements between each of the masses (Fig. 21). We provide details of a procedure for obtaining the dispersion curves as described in Refs. [340–342]. Considering periodicity, the set of homogeneous equations describing the motion of each mass in the unit cell is obtained as follows:

Fig. 21
Fig. 21
Close modal
$m1u··1j+(c1+c2)u·1j-c2u·2j-c1u·2j-1+(k1+k2)u1j -k2u2j-k1u2j-1=0$
(54)
$m2u··2j+(c1+c2)u·2j-c2u·1j-c1u·1j+1+(k1+k2)u2j-k2u1j-k1u1j+1=0$
(55)

where $uαj$ is the displacement of mass α in an arbitrary jth unit cell. In this notation, a unit cell and its neighbors may be identified by j + r, where r = 0, −1, +1 denote the present, previous, and subsequent unit cells, respectively.

As explained earlier, it is necessary to adopt a generalized form of Bloch's theorem—in which the frequencies are allowed to be complex—in order to obtain a dispersion relation that represents the free wave propagation characteristics of the chain. For the 1D discrete model, the generalized Bloch wave solution takes the form
$uαj+n(x,κ;t)=U∧αeλt=U¯αei(κx+nκa)=U˜αei(κx+nκa)+λt$
(56)
where $U∧α,U¯α$, and $U˜α$ represent different forms of the complex wave amplitude and λ a complex frequency. Applying Eq. (56) to Eqs. (54) and (55), and arranging in matrix form, yield the following complex eigenvalue problem:
$[A11A12A21A22][U˜1U˜2]=[00]$
(57)
where
$A11=λ2m1+λ(c1+c2)+k1+k2,A12=-λ(c1e-iκa+c2)-(e-iκak1+k2),A21=-λ(c1eiκa+c2)-(eiκak1+k2),A22=λ2m2+λ(c1+c2)+k1+k2$
which may also be expressed as
$[λ2M+λC(κ)+K(κ)]U˜=0$
(58)
where $U˜=[U˜1 U˜2]T$, and the mass matrix M, damping matrix C(κ), and stiffness matrix K(κ) are, respectively, defined as
$M=[m100m2],C(κ)=[c1+c2-(c1e-iκa+c2)-(c1eiκa+c2)c1+c2],K(κ)=[k1+k2-(k1e-iκa+k2)-(k1eiκa+k2)k1+k2]$
In order to solve Eq. (58) for the roots, we recast the problem into a first-order state-space formulation by the transformation of variables, $Y¯=[¯·U U¯]T$, to get
$[0MMC(κ)]¯·Y+[-M00K(κ)]Y¯=0$
(59)
The Bloch solution, in turn, is also written in terms of the transformed variables, that is,
$Y¯=Y˜eγt$
(60)
Upon substituting Eq. (60) into Eq. (59), we get
$([0MMC(κ)]γ+[-M00K(κ)]Y˜)=0$
(61)
which is now a standard complex eigenvalue problem. Given their orthogonality with respect to the system matrices, the eigenvectors decouple the equations into four modal equations with complex roots $γs*,s*=1,…,4$ appearing in complex conjugate pairs and thus effectively representing two (s = 1, 2) single degree-of-freedom systems. Thus, we can write
$γs=-ζsωs±iωds=-ζsωs±iωs1-ζs2, s=1,2$
(62)
If we focus our attention on only the first eigenvalue in each complex conjugate pair, then we can extract the damped wave frequency as
$ωds(κ)=Im[γs(κ)], s=1,2$
(63)
and the corresponding damping ratio
$ζs(κ)=-Re[γs(κ)]Abs[γs(κ)], s=1,2$
(64)

Note that $ωs(κ)=Abs[γs(κ)]$. For the special case of proportional (Rayleigh) viscous damping, the value of ωs is equal to the undamped frequency [343].

To get some insight into the effects of damping on the wave propagation characteristics of phononic materials and structures, we generate dispersion curves for the 1D chain as shown in Fig. 21. In addition to customary frequency dispersion curves—which for our damped model we obtain from Eq. (63)—we may also plot curves for the Bloch mode damping ratios using Eq. (64). For demonstration, a specific set of material parameters are used: rm = m2/m1 = 9, rc = c2/c1 = 0.5, rk = k2/k1 = 1, and $ω¯=k2/m2=149.07 rad/s$. A damping parameter, β = c2/m2, is varied to show the dependence on the damping intensity. In Fig. 22, a frequency band diagram appears on the top and a damping ratio band diagram appears on the bottom. While these results are specific to the values chosen for the parameters rm, rc, and rk, they provide a basic insight on the effects of damping on the elastodynamic behavior of a periodic material or structure.

Fig. 22
Fig. 22
Close modal

We observe in Fig. 22 (top) that as the damping intensity $β/ω¯$ is increased, the optical branch drops while the acoustic branch experiences little change—this in turn leads to a reduction in the size of the bandgap. It also leads to a branch-overtaking phenomenon [340–342] whereby a higher branch overtakes a lower branch as damping is raised. Furthermore, we observe that the descent of the optical branch takes place at slightly faster rates at low wavenumbers compared with high wavenumbers. The damping ratio diagram in Fig. 22 (bottom) follows a consistent trend with an indication that the effect of damping is slightly more significant at low wavenumbers. In looking further at both diagrams in tandem, we observe that the optical branch cuts-off for long wavelengths that are overdamped, and the cut-off wavenumber is precisely at ζ2 = 1 (which corresponds to the critical damping ratio) [341,342]. This cut-off phenomenon implies that under high damping conditions, there is a possibility for wavenumber bandgaps to arise as illustrated in Fig. 23, which displays a bandgap map for both frequencies and wavenumbers.

Fig. 23
Fig. 23
Close modal
Fig. 24
Fig. 24
Close modal

#### Damped 1D Lattice With Internal Resonator.

In this section, we present the effects of damping on the band structure of a 1D locally resonant diatomic chain, as depicted in Fig. 24. Following the same analysis procedure described in the Sec. 2, we derive the equations of motion for this system leading to the matrix equation given in Eq. (58), except now the matrices are defined as follows:
$M=[m100m2],C(κ)=[2c1(1-cosκa)+c2-c2-c2c2],K(κ)=[2k1(1-cosκa)+k2-k2-k2k2]$
The frequency dispersion curves and damping ratio curves for this 1D locally resonant lattice (using the same mass, stiffness, and damping properties as in Sec. 4.1.2) are depicted in Fig. 25,
Fig. 25

Frequency (top) and damping ratio (bottom) band structures for the 1D viscously damped locally resonant lattice depicted in Fig. 24. Acoustic branches are in solid lines and optical branches are in dashed lines (adapted from Ref. [342] with permission from Springer).

Fig. 25

Frequency (top) and damping ratio (bottom) band structures for the 1D viscously damped locally resonant lattice depicted in Fig. 24. Acoustic branches are in solid lines and optical branches are in dashed lines (adapted from Ref. [342] with permission from Springer).

Close modal
. Contrasting these results with those observed in Fig. 22 highlight the qualitative difference between the effects of damping on a phononic crystal and a locally resonant metamaterial.

The most striking contrast is that in the latter case, the drop rate of the optical branch is significant at high wavenumbers compared with low wavenumbers. This is reflected in the rapid rise of the dissipation factors at high wavenumbers as shown at the bottom of Fig. 25. We also observe that at low wavenumbers, the damping ratio values corresponding to the optical modes for the conventional lattice exceed those of the locally resonant lattice for a given damping intensity $β/ω¯$ and vice versa at high wavenumbers. With regard to the damping ratios of the acoustic branch modes, these are similar at low wavenumbers but are higher for the regular lattice at high wavenumbers.

The qualitative and quantitative differences revealed in the dissipation diagrams as shown at the bottom of Figs. 24 and 25, respectively, have recently been examined in more detail [344]. In order to yield some insights into the difference in the overall damping capacity, a comparison was made between the same two models except with a slight alteration in the parameter values in order to allow both systems to exhibit an identical value of long-wave speed. It was shown that despite both systems having the same values of prescribed damping in the dashpots, and being statically equivalent, the locally resonant lattice (i.e., the metamaterial) exhibits more dissipation across the spectrum. This phenomenon, which has been termed metadamping, is indicative of an emergence of damping due to the presence of local resonance and may be utilized toward the design of materials that need to be simultaneously stiff and dissipative.

#### Damped 2D Phononic Crystal.

We now turn to the treatment of damping in the study of Bloch wave propagation in multidimensional continuous media, such as a phononic crystal. The starting point for tackling this problem is the standard balance law for elastodynamic forces
$∇·σ=ρu··$
(65)
where $σ$ is the stress tensor. For the constitutive behavior, we will assume linear and isotropic elastic response and some general type of damping
$σ=C:∇su+αd(u·)$
(66)

where C is the conventional elasticity tensor, $αd$ is the velocity-dependent dissipative component of the stress tensor.

Substituting Eqs. (66) and (35) into Eq. (65) yields
$∇·C:∇su+αd(u·)=ρu··$
(67)
At this point, the generalized form of Bloch's theorem is utilized, as done for the discrete case
$u(x,κ;t)=u¯(x;t)ei(κTx), u¯(x;t)=u˜(x)eλt$
(68)
where $u¯$ and $u˜$ are the time-dependent and time-independent displacement Bloch functions, respectively. We emphasize that the generalized form of Eq. (68) allows for both spatial and temporal attenuation in the continuous Bloch wavefield. The choice of Eq. (68) for the solution has been validated in Refs. [340,345,346] by conducting a comparison between the location of bandgaps exhibited by a damped phononic crystal and corresponding natural frequency gaps for a finite structure composed of several unit cells of precisely the same construction as the phononic crystal. The results have shown that the two sets of gaps match perfectly. The choice of damping model for $αd(u·)$ is abundant. For simplicity, we will choose viscous damping—the same type of damping we adopted in Secs. 4.1.2 and 4.1.3. At this point, it is practical to discretize the problem using a numerical technique, such as the finite element method (see Sec. 3.1.3 for details). Upon substitution of Eq. (68) into Eq. (67), application of periodic boundary conditions, and discretization into an N degree-of-freedom system, we obtain
$[λ2M+λC(κ)+K(κ)]U˜=0$
(69)
where M, C, and K now denote the global finite element mass, damping, and stiffness matrices, respectively. Inspection of Eq. (69) immediately reveals a resemblance to Eq. (58) with the exception that the independent variable is now $κ$, which is a wavevector rather than a wavenumber. This congruence of format allows us to mechanically apply the steps described in Sec. 4.1.2 following Eq. (58) to the present problem considering that only the size and components of the system matrices differ. Finally, we obtain the eigensolution that embodies both the frequency of damped wave motion ωds and the damping ratio ζs for mode s
$λs(κ)=-ζs(κ)ωs(κ) ± iωds(κ), s=1,2,…,N$
(70)
We now examine the behavior of a continuous model of a generally damped phononic crystal. While the formulation is applicable to three-dimensional models, we limit our attention in this example to a 2D plain-strain model. Specifically, a square lattice is considered with a bimaterial unit cell consisting of a centrally located square inclusion (with a filling ratio of 0.3086), as shown in the inset of Figs. 26,
Fig. 26

Effect of viscous damping on wave dispersion (case 1): (a) frequency dispersion diagram and (b) damping ratio diagram. In both subfigures, the undamped and damped models are represented by black lines and red dashed lines, respectively.

Fig. 26

Effect of viscous damping on wave dispersion (case 1): (a) frequency dispersion diagram and (b) damping ratio diagram. In both subfigures, the undamped and damped models are represented by black lines and red dashed lines, respectively.

Close modal
and 27,
Fig. 27

Effect of viscous damping on wave propagation directionality: Isofrequency contour lines for model without (black lines) and with (red lines) damping. The topology of the unit cell considered is shown in the inset (white: compliant/low density material; black: stiff/high density material).

Fig. 27

Effect of viscous damping on wave propagation directionality: Isofrequency contour lines for model without (black lines) and with (red lines) damping. The topology of the unit cell considered is shown in the inset (white: compliant/low density material; black: stiff/high density material).

Close modal
. The material phase for the matrix (denoted by subscript “1”) is chosen to be compliant and light, while the phase for the inclusion (denoted by subscript “2”) is stiff and dense. The ratio of the Young's moduli is E2/E1 = 20 and that of densities is ρ2/ρ1 = 2, and a Poisson ratio of 0.34 is assumed for both phases. As for damping, we consider for simplicity a proportionally damped model, where $C(κ)=pM+qK(κ)$, and examine two cases. In case 1, we take p = 0 and q = c*/10a and in case 2, we take p = 2a/c* and q = c*/10a ($c*=E1/ρ1$ and a denotes the unit-cell side length, or lattice constant, as in the discrete problems covered earlier). We use the finite element method to discretize the spatial domain into a mesh consisting of 9 × 9 uniformly sized four-node bilinear quadrilateral elements. As described above, we obtain the eigensolution and extract $ζs(κ)$ and $ωds(κ)$ from Eq. (70). Figure 26(a) shows the frequency dispersion curves for the damped phononic crystal (case 1) and, for comparison, a matching undamped phononic crystal for which p = q = 0. The corresponding damping ratio curves are shown in Fig. 26(b). Figure 27 shows isofrequency contours (corresponding to the second branch) as a function of the x- and y-components of the wavevector $κ$ for case 2. The results demonstrate the remarkable effect that damping has on the frequency and directionality of wave propagation. This effect may be utilized in designing waveguides using damped phononic materials or structures.

While the state-space approach presented in this section provides an exact treatment of damping (unless a numerical discretization technique is applied), it often faces numerical conditioning issues due to the loss of symmetry in the first-order eigenvalue problem. This is more likely to occur for relatively high values of prescribed damping (i.e., high enough for a particular branch to approach the critical damping level). The recently proposed Bloch–Rayleigh perturbation method may be used to overcome this limitation [347]. This method utilizes the Bloch waves of an undamped version of the phononic medium as basis functions to provide approximate closed-form expressions for the complex eigenvalues and eigenvectors of the damped system.

### Nonlinear Phononic Materials and Structures.

In this section, we first review studies concerned with nonlinear response in phononic materials and structures before providing a detailed technical review of a perturbation approach for predicting dispersion behavior in weakly nonlinear periodic systems. Also highlighted is the potential for nonlinear periodic systems to enable tunable filters, waveguides, resonators, frequency isolators, and acoustic diodes.

#### Literature Review.

The number of investigations carried out on nonlinear phononic materials and structures is quite limited in comparison with linear investigations. The influence of nonlinearities on the propagation and attenuation zones in chain-like one-dimensional periodic systems with weak nonlinearities was the focus of early studies. Vakakis and his co-authors [348] analyzed nonlinear chains subjected to external forcing and ground springs using a long wavelength continuum approximation together with the multiple scales perturbation method. They identified propagation and attenuation zones and discussed localized response. Next, Vakakis and King [349] studied a chain composed of monocoupled continuous subsystems with cubic nonlinearities using the multiple scales approach and presented the linear band diagrams before discussing synchronous or standing wave solutions in the attenuation zones. Their treatment of propagation zones was much shorter than their attenuation zone treatment; however, they did develop nonlinear dispersion relationships that exhibit amplitude dependence. Chakraborty and Mallik [350] studied the cubic chain and assumed a Bloch solution at the outset of their analysis and then studied the effect of nonlinearities on the propagation constant and natural frequencies. They used their approach to find nonlinear normal modes. Lazarov and Jensen [351] considered a linear chain with attached nonlinear oscillators. Their analysis employed a harmonic balance approach to uncover amplitude-dependent dispersion. Wang and Wang recently considered a similar problem where the linear chain was replaced by a continuous beam [352]. Another recent investigation involving the combination of nonlinearity and local resonance utilized a closed-form finite-strain dispersion relation in conjunction with the transfer matrix method [353]. Marathe and Chatterjee [354] looked at a damping nonlinearity and used harmonic balance and multiple scales to uncover the decay rate in the propagation zone. Narisetti et al. [355] developed a Lindstedt–Poincaré perturbation technique, valid for all wavelengths, for analyzing cubic chains with an arbitrary number of masses in their unit cells. The specific intent of their work was to capture dispersion and bandgap shifts, and to then suggest novel devices this might enable, to include frequency isolators and tunable switches and filters. Manktelow et al. [356] studied similar systems, this time with a multiple scales approach, and documented wave–wave interactions and a new dispersion shifting term due to the presence of multiple waves at different frequencies/wavenumbers.

Strongly nonlinear 1D systems have also been studied from the standpoint of phononic systems. Nesterenko [357] used a continuum approximation to derive exact solutions to chain-like granular media with interactions governed by a Hertzian contact law. In this chain, the stiffness goes to zero when masses separate, making it very difficult to use asymptotic approaches. Nesterenko documented soliton behavior characterized by amplitude-dependent phase velocity and reported on collisions of multiple solitons. Building upon Nesterenko's work, Chatterjee [358] developed a new asymptotic solution for the traveling wave that is valid near the wave crest. Other authors [359–362] report experimental results confirming the existence of solitary waves in granular media, and also report tunability of wave speeds and bandgaps with precompression. Spadoni and Daraio [164] propose using such systems for generation and control of “sound bullets.” Boechler et al. [363] studied discrete breathers (i.e., localized solutions at frequencies in the attenuation zone) in a 1D diatomic granular crystal. Starosvetsky and Vakakis [364] presented an analytical study of an uncompressed chain of beads, without invoking the continuum limit, and found families of stable multihump traveling-wave solutions with arbitrary wavelengths. They also found strongly localized out-of-phase standing waves. Narisetti et al. [365] focused on plane waves in strongly nonlinear periodic media and analyzed their propagation using the harmonic balance method. Liang et al. [165] described a novel acoustic diode device from the coupling of a superlattice with a strongly nonlinear medium.

Two-dimensional nonlinear phononic systems have also been studied. Discrete breathers in 2D lattices were studied by Franchini et al. [366] and Feng and Kawahara [367]. Soliton analysis in 2D lattices date back to the seminal study by Toda [368]. Soliton studies since Toda's paper have been numerous. A continuum treatment of soliton behavior was presented by Sreelatha and Joseph [369]. Narisetti et al. [370] generalized their Lindstedt–Poincaré approach to any dimensionality and reported amplitude-dependent group velocities and tunable wave beaming. Prestress as a means of achieving tunability has also been studied. For example, Gilles and Coste [371] experimentally studied low and high-frequency behavior of wave propagation in a triangular lattice of spherical beads under isotropic stress. Bertoldi, Boyce, and their co-authors studied tuning of phononic properties via prestress [192] and opening of new bandgaps [372].

#### Band Structure Calculation by Perturbation Analysis.

In this section, we provide a technical review of a recently presented perturbation approach for analyzing wave propagation in weakly nonlinear phononic systems [355,370]. The approach is valid for discrete systems propagating both long and short wavelength waves, i.e., it is not necessary to invoke a continuum approximation. For analysis of continuous systems, the approach is applied following a spatial discretization of the problem (e.g., lumped parameter, finite element, etc.). As in a linear analysis, the goal is to determine the system's band structure by analyzing a unit cell, thus reducing the infinite system to a finite one. This is complicated by the presence of nonlinear terms in the equations of motion, which do not allow Bloch's theorem to be invoked in the usual manner [373]. This is overcome in the perturbation analysis by asymptotically expanding the equation set and moving nonlinear terms to higher orders where they can be treated as inhomogeneous forcing terms. Thus, the governing equations at each order have linear kernels amenable to Bloch analysis.

We start by considering periodic 1D, 2D, and 3D phononic systems that can be in either lumped or continuous form. The specific case of a continuous system discretized using the finite element method is considered since it provides all of the generality necessary to analyze almost any weakly nonlinear phononic system, and since finite element-specific issues arise in the procedure that must be considered. The systems of interest are depicted in Fig. 28 and consist of a central unit cell (colored red) surrounded by nearest-neighbor unit cells, which are in turn acted upon by external forces f arising from next-nearest neighboring unit cells. This specific assembly of unit cells is chosen such that the central unit cell does not experience external forces. Following introduction of the perturbation technique, the nearest-neighbor unit cells are removed from consideration using appropriate application of Bloch's theorem.

Fig. 28
Fig. 28
Close modal
Suppose that the continuous domain has been discretized using a finite-element technique such that qi denotes the DOFs for the ith node. Furthermore, let $u(p,q,r)=[q1,q2,…qn]T$ denote the collection of DOFs within the unit cell indexed by p, q, and r (only p for a 1D system, or only p and q for a 2D system). We choose to order the collection of generalized coordinates for a given unit cell and its nearest neighbors as
$usys=[u(-1,-1,-1),u(-1,-1,0),…,u(0,0,0),…,u(0,1,1),u(1,1,1)]T$
(71)
Given this ordering scheme, the discretized equations of motion can be written as
$Msysu··sys+Ksysusys=fsys$
(72)
where Msys and Ksys denote the mass and stiffness matrices for the central unit cell and its neighbors as shown in Fig. 28, and fsys denotes both the external forces acting on the nearest-neighbor unit cells and the nonlinear contribution of internal forces throughout the assembly. The resulting dynamic equations constitute an open set of nonlinear difference equations whose solution is dependent upon all u(p,q,r) surrounding the central unit cell. By lumping the mass matrix [374] for each unit cell, the equations of motion governing the central unit cell can be isolated from the others and are given by
$Mu··(0,0,0)+∑p=-11∑q=-11∑r=-11K(p,q,r)u(p,q,r)=f(0,0,0)$
(73)
where M denotes a partition of the total mass matrix corresponding to the central unit cell, K(p,q,r) denote partitions of Ksys, and f(p,q,r) denotes the partition of the system force vector holding only nonlinear terms (external forces on the central unit cell are identically zero by choice of assembly and absence of external loading). Note that the isolation procedure above can be visualized as retaining only the rows of equations containing the central unit cell's mass matrix. Assuming weakly nonlinear behavior, the force vector f(p,q,r) can be scaled
$f(p,q,r)=ɛf(p,q,r)NL$
where $|ɛ|≪1$ is a small parameter quantifying the magnitude of the nonlinearity. For a general 3D system, the open set of difference equations (Eq. (73)) for the central unit cell then appears as
$Mu··+∑p=-11∑q=-11∑r=-11K(p,q,r)u(p,q,r)=ɛfNL$
(74)

where the simplified notation u(0,0,0) = u and $f(0,0,0)NL=fNL$ has been employed.

The small parameter ε appearing in Eq. (74) facilitates an asymptotic solution approach. Herein, we review an approach based on the Lindstedt–Poincaré method, whereas an approach based on multiple scales can be found in Ref. [356]. We begin by introducing a dimensionless time τ = ωt and standard expansions
$ω=ω0+ɛω1+O(ɛ2)u(p,q,r)=u(p,q,r)(0)+ɛu(p,q,r)(1)+O(ɛ2)$
(75)
These expansions, when substituted into Eq. (74), yield ordered equations
$O(ɛ0): ω02M∂2u(0)∂τ2+∑p=-11∑q=-11∑r=-11K(p,q,r)u(p,q,r)(0)=0$
(76)
$O(ɛ1): ω02M∂2u(1)∂τ2+∑p=-11∑q=-11∑r=-11K(p,q,r)u(p,q,r)(1) =-2ω0ω1M∂2u(0)∂τ2+fNL(u(p,q,r)(0))$
(77)
Since the discretized system at O2) is periodic and linear, Bloch's theorem may be applied to Eq. (76) to produce a set of equations describing an infinite linear system. This equation set depends on the dimensionless Bloch wavevector $μ=[μ1,μ2,μ3]$
$ω02M∂2u(0)∂τ2+K(μ)u(0)=0$
(78)
where the wavenumber-reduced stiffness matrix $K(μ)$ is
$K(μ)=∑p=-11∑q=-11∑r=-11K(p,q,r)e-i(pμ1+qμ2+rμ3)$
Seeking solutions that are time-harmonic with amplitude A
$u(0)=A2φjeiτ+c.c.$
with c.c. denoting a complex conjugate, Eq. (78) yields a standard eigenvalue problem for eigenvalues $ω0,j2(μ)$ and Bloch wave modes $φj(μ)$, parameterized by the dimensionless Bloch wave vector
$(K(μ)-ω02M)φj=0$

The dispersion relation for the continuous system is approximated by $ω0,j(μ)$, where the jth eigenvalue corresponds to the jth dispersion branch (or surface). Note that the approximate dispersion relations become more accurate as the number of mesh elements is increased.

Expanding the nonlinear vector $fNL(u(0))$ from Eq. (77) in a Fourier series
$fNL=∑n=-∞n=∞cj(A)einτ$
(79)
and keeping only the c1 terms, the O1) system of equations (77) may be reduced by assuming a “secular-Bloch” type solution (see Ref. [370] for details). The eigenvectors, or Bloch wave modes, for the reduced O0) problem are the same as those for the reduced O1) homogeneous system, and thus the matrix of wave modes $Φ=[φ1...φn]$ decouples the system of equations in a similar manner. It can be shown that removal of secular terms leads to an O1) frequency correction term for the jth dispersion branch (or surface) given by
$ω1,j(μ)=-φjHc1ω0,jAφjHMφj$
(80)
The updated dispersion relation is then given as
$ωj(μ)=ω0,j(μ)+ɛω1,j(μ)+O(ɛ2)$

which may be used to identify amplitude-depending shifting of dispersion bandgaps. It is clear that the Bloch wavemodes $φj(μ)$ and frequencies $ω0,j(μ)$ in Eq. (80) are readily obtained from the mass matrix M and wavenumber-reduced stiffness matrix $K(μ)$ of the O0) problem. However, the first Fourier coefficient c1 requires additional consideration.

Simple systems, such as the nonlinear diatomic spring-mass system, admit tractable analytical expressions for fNL and thus c1 [355]. However, for more complex systems, such as those with inclusions, voids, nontrivial geometry, or complex nonlinearities, hand calculations for explicit expressions become impractical. Instead, numerical estimates of c1 can be obtained using commercial software such that the complex design space of periodic materials can be investigated. More details on 1D implementation can be found in Ref. [375].

#### Nonlinear Behavior in Example Systems.

The nonlinear asymptotic approach is applied next to an example one-dimensional system and illustrative behavior is discussed. This example considers a one-dimensional diatomic chain with a cubic stiffness nonlinearity, which exhibits amplitude-dependent dispersion shifts and bandgaps. This behavior can be exploited, for example, to devise a nonlinear frequency isolator, as described at the end of the section. Two- and three-dimensional nonlinear systems offer additional device implications, arising from amplitude-dependent group velocities, but are not discussed herein.

Linear diatomic chains behave as stop band filters in which a center frequency band does not permit wave propagation. Note that as the plane-wave frequency approaches the acoustic cut-off frequency, the lighter masses in the chain rest. Alternatively, as the frequency approaches the optical cut-off frequency, the heavier masses rest. Due to this behavior and the need to define ratios of mass amplitudes, the perturbation approach must be applied to each branch separately to avoid nonuniformity in the asymptotic expansions over all frequencies. Therefore, frequency corrections for each mode are expressed in terms of the amplitude of the masses that do not rest as the excitation frequencies tend toward the cut-off zone. This is discussed again below when introducing ηopt and ηaco.

The diatomic chain assembled mass, stiffness, and nonlinear force vectors required in Eq. (74) are given by
$M=[m100m2], K(-1)=[0-k00],K(0)=[2k-k-k2k], K(1)=[00-k0],fNL=-[Γ(u(0),1-u(-1),2)3+ Γ(u(0),1-u(0),2)3Γ(u(0),2-u(1),1)3+ Γ(u(0),2-u(0),1)3]$
where u(p),l denotes the lth degree of freedom in the pth unit cell, and further it is noted that K(p,q,r) reduces to K(p) since the problem is one-dimensional. The wavenumber reduced stiffness matrix then follows as
$K(μ)=[2k-k(1+eiμ)-k(1+e-iμ)2k]$
At O(ε0), the linear dispersion relationships are recovered with the acoustic and optical eigenvalues given by
$ω0aco=ωn1+β-1+β2+2βcos(μ),ω0opt=ωn1+β+1+β2+2βcos(μ)$
where β ≡ m1/m2 is assumed less than one without loss in generality, $ωn≡k/m1$, and numeral subscripts are replaced by superscript aco and opt to more easily convey physical significance. The eigenvectors can be expressed as the ratio of one element to the other
$ηaco≡φ1acoφ2aco=1+eiμ1-β+1+β2+2βcos(μ),ηopt≡φ2optφ1opt=1-β-1+β2+2βcos(μ)1+eiμ$

where $φkaco$ denotes the kth element in $φaco$; similarly for $φkopt$.

The remaining expressions for c1 and the corrections $ω1aco$ and $ω1opt$ are difficult to express in a compact form. However, results from the perturbation procedure are presented in Fig. 29 for an example parameter set corresponding to a hardening nonlinearity. Note that the frequency corrections shift both dispersion branches upwards with increasing amplitude—a softening nonlinearity would do the opposite. The corrections preserve the property of zero group velocity at the edge of the Brillouin zone (μ = π) for both branches, as well as the zero group velocity property of the optical branch at μ = 0. Finally, the corrections can be observed to increase the center frequency and width of the bandgap in this example.

Fig. 29
Fig. 29
Close modal

The amplitude-dependent shifts observed in systems, such as the nonlinear diatomic chain can be exploited in new ways to create tunable devices. These can include tunable resonators, filters, and waveguides that operate under the action of signal gain and use amplitude-dependent bandgaps to achieve tunability. They can also include devices that manipulate frequencies in novel ways. As an example of the latter, a nonlinear frequency isolator (see Fig. 30) was proposed in Ref. [355], which uses a tunable signal gain to control whether one or two frequencies pass through the isolator. The input signal contains two frequencies, one at 1.75 rad/s lying significantly inside of the acoustic propagation zone, and the other at 2.46 rad/s lying in the optical propagation zone at low amplitudes and outside of it at high amplitudes; see Fig. 30(b). Numerical simulations confirm the operation of the device. Figure 30(c) demonstrates that at low gain corresponding to a wave amplitude A = 0.5, both frequencies propagate and pass through the device. Figure 30(d) demonstrates that at higher gain, corresponding to a wave amplitude A = 0.75, only the lower frequency passes through the device.

Fig. 30
Fig. 30
Close modal

Two-dimensional nonlinear systems demonstrate amplitude-dependent spatial characteristics which that further increase the design space of nonlinear phononic systems. In addition to amplitude-dependent bandgap behavior seen in 1D systems, group velocity, and hence energy propagation, also exhibits amplitude dependence. This can be used, for example, to design spatial filters, multiplexors, and other logic functions based on amplitude-dependent wave steering. Interested readers can consult Refs. [365,370,375] for further details. In all these applications, as in the linear problem, the unit-cell topology could be optimized for desired performance [376].

## Experimental Studies

This last section is devoted to an overview of examples of experimental activities conducted in the field of phononic materials and structures. Following this brief review of the literature, most of the section is devoted to the presentation of experimental configurations investigated by one of the authors. The examples presented illustrate simple measurement techniques applied for the analysis of the considered class of systems and document physical evidence of phenomena discussed in Secs. 2–4.

The number of papers reporting experimental results on the dynamics of periodic structures is significantly lower than those dealing with theoretical and numerical aspects. However, recently there has been an increase in the number of publications with experimental results illustrating phenomena, such as bandgaps, negative refraction, attenuation due to internal resonances, and directional acoustic radiation. This observation particularly applies to the phononic crystals and metamaterials literature, where arrays of metallic cylinder submerged in water, stubbed plates, and micromachined lattice structures are often considered as test samples [214,377–383].

In structural dynamics, early studies have investigated the modal properties of stiffened plates and shells as examples of structural components found in many engineering constructions. In particular, in Ref. [384], attention is devoted to relating the modal density of the structural assembly with the resonances of the unit cell, or single bay, and to the mapping of the theoretical stop and pass bands onto measured spectra and sonograms. Of immediate practical relevance is the investigation of the noise radiation characteristics of underwater stiffened shells as reported, for example, in Refs. [385–390]. In particular, the scattering spectrum of ribbed shells is explored in Ref. [385], where confirmation is provided that the flexural Bloch waves, and Bragg scattering from the ribs, are principle sources of scattering for typical ribbed shells in a frequency range related to the ring frequency. The acoustic scattering of shells with multiple internal resonators has been investigated in Refs. [387,388], where the analysis in the frequency/wavenumber domain provides an effective visualization of the band structure [390]. Other examples of experimental studies in structural dynamics include the investigation of a vibration isolation scheme for beams featuring a periodic structure [391], and the characterization of directional response in a beam grillage [77].

In one of the early studies on sonic or phononic materials, the case of a layered structure consisting of water and a perspex plate is studied theoretically and experimentally [392]. In particular, the experimental results show clear evidence of bandgaps, along with the effect of a vacancy defect which produces a narrow passband in each of the bandgaps. Experimental investigations of bandgaps in micromachined air/silicon phononic structures are presented in Ref. [104], while Ref. [393] illustrates experimentally (and theoretically) that super resolution can be achieved while imaging with a flat lens consisting of a phononic crystal exhibiting negative refraction. Experimental demonstration of negative refraction of transverse elastic waves is reported for a 2D phononic crystal made of a square lattice of cylindrical air cavities in an aluminum matrix in Ref. [301]. Similar demonstrations are presented in Ref. [383], while evidence of bandgaps in 2D systems is presented in Refs. [377,378,381]. In the GHz range, experimental ripple wavefield visualization in a phononic crystal is reported by Profunser and coworkers [394,395]. In conclusion to this short, and nonexhaustive, review of experimental work, it is worth mentioning the investigation of localization in addition to bandgaps in one-dimensional systems as illustrated in Ref. [382].

### Stubbed Plate for Low-Frequency Wave Attenuation.

A first example illustrates low-frequency attenuation achieved through the resonant behavior of oscillators connected to the primary structure. The concept is illustrated in Sec. 2.3 and has been exploited for low-frequency/subwavelength attenuation as reported in Ref. [396], or to achieve negative mass effects, as discussed in Refs. [224,240,243,296].

Here, we illustrate the experimental investigation of a plate with mechanical resonators tuned at frequencies ranging from 650 to 3500 Hz. These frequencies are significantly lower than those corresponding to Bragg scattering based on the spatial period of the array that is of the order of 100 kHz. In the considered configuration, originally presented in Ref. [397], a thin aluminum plate features a periodic array of composite stubs made of tungsten and silicone rubber (Fig. 31). The spatial period of the array and the diameter of the stubs are selected, respectively, as a = 1 cm and d = 6 mm, while the stubs are h = 3 mm high. The plate is excited by a surface mounted PZT transducer, and the resulting wavefield is recorded by a scanning laser Doppler vibrometer (SLDV) over a 2D measurement grid.

Fig. 31
Fig. 31
Close modal
The SLDV data corresponding to the displacements along the directions of symmetry of the structure (horizontal, vertical, diagonal), denoted as $u∧θ(r,f)$ (with θ = 0 deg, 45 deg, 90 deg indicating the considered direction shown in Fig. 31), are extracted from the recorded field, $u∧(x1j,x2j,f)$ (with x1j, x2j denoting the coordinates of the jth point of the measurement grid), according to the following averaging procedure:
$u∧r(r,f)=1N∑j|(x1j,x2j)∈θu∧(x1j,x2j,f)$
(81)
Transmission data are evaluated by normalizing the plate response to the displacement recorded at the center of the array, i.e., $u∧0(f)=u∧(0,0,f)$. Results are presented in Fig. 32,
Fig. 32

Transmission measurements on plate with resonant stubs show strong evidence of bandgap behavior in the 1300–3500 Hz (Reprinted with permission from Ref. [397]. Copyright 2012, AIP Publishing LLC).

Fig. 32

Transmission measurements on plate with resonant stubs show strong evidence of bandgap behavior in the 1300–3500 Hz (Reprinted with permission from Ref. [397]. Copyright 2012, AIP Publishing LLC).

Close modal
, where strong attenuation is observed in the 1300–3500 Hz range. The behavior of the plate is further illustrated by presenting plots of its dynamic deformed shapes at selected frequencies. These plots are obtained by taking frequency snapshots of the recorded wavefield from its frequency-domain representation. Figure 33,
Fig. 33

Measured steady-state deformed shapes, representing the stationary response of the plate for excitation at 0.79 kHz (below the bandgap) (a), 2.35 kHz (inside the bandgap) (b), and 5.39 kHz (above the bandgap) (c). The red circle represents the source excitation of the acoustic waves (Reprinted with permission from Ref. [397]. Copyright 2012, AIP Publishing LLC).

Fig. 33

Measured steady-state deformed shapes, representing the stationary response of the plate for excitation at 0.79 kHz (below the bandgap) (a), 2.35 kHz (inside the bandgap) (b), and 5.39 kHz (above the bandgap) (c). The red circle represents the source excitation of the acoustic waves (Reprinted with permission from Ref. [397]. Copyright 2012, AIP Publishing LLC).

Close modal
shows three dynamic deformed shapes corresponding to frequencies below the bandgap (Fig. 33(a)), inside the bandgap (Fig. 33(b)), and above the bandgap (Fig. 33(c)). At the frequency inside the bandgap, the plate appears virtually undeformed because the input energy is absorbed by the composite resonators and remains confined near the source. However, for excitation below and above the bandgap, the plate deforms globally according to its modal response at the considered frequency.

Waveguiding is also investigated by creating a linear defect by removing three rows of stubs, which corresponds to a waveguide of about 2.4 cm width. The results of Fig. 34 clearly illustrate the confinement of the wavefield within the designed waveguide, for excitation within the bandgap. The excellent energy confinement obtained through the linear defect is very evident, and it is particularly important to observe how the wavelength of deformation is noticeably larger than the spatial spacing of the stubs, to indicate that such waveguiding is not the result of Bragg scattering.

Fig. 34
Fig. 34
Close modal

### Beam With Periodic Electromechanical Resonators.

In the next example, we illustrate a 1D metamaterial with internally resonating units consisting of periodic arrays of piezoelectric patches bonded to a 1D waveguide, i.e., a beam in bending (see Fig. 35). A detailed description of the experiments can be found in Ref. [289]. The piezoelectric patches are shunted through a resonating, passive circuit, and therefore act as equivalent resonators of the kind illustrated in Sec. 2.3.

Fig. 35
Fig. 35
Close modal
According to Ref. [398], shunting of the piezoelectric patch with electrodes across its thickness direction modifies the Young's modulus of the shunted patch based on the following expression:
$EpSU(ω)=EpD(1-k3121+iωCpɛZSU(ω))$
(82)

where $Cpɛ$ the electrical capacitance of the piezo at constant strain and ZSU(ω) is the electrical impedance of the shunting circuit. Also in Eq. (82), k31 denotes the electromechanical coupling coefficient and $EpD$ is the Young's modulus of the piezoelectric material when the shunting network is in an open circuit configuration ($ZSU(ω)→∞$). This quantity is related to the piezo Young's modulus Ep through the expression $EpD=Ep/(1-k312)$. When a resistor-inductor shunt is applied to the piezoelectric material, $ZSU(ω)=R+iωL$, the piezo Young's modulus features a resonant behavior at a frequency ωT that can be selected by tuning the inductor $L=1/(ωT2Cpɛ)$ [398]. Therefore, such resonant frequencies can be conveniently tuned through proper selection of the electrical impedance connected to each patch, so that no modification to the structure is necessary.

The dispersion characteristics, bandgaps, and wave attenuation properties of the beam are evaluated experimentally through wavefield measurements of the kind presented in Sec. 5.1. The recorded data have fine temporal and spatial resolution, which is sufficient for the estimation of the frequency/wavenumber content of the beam response through the evaluation of the spatial/temporal 2D FT.

Tests were conducted on the beam shown in Fig. 36. The beam, which is made of aluminum and is 1.6 m long, features an array of 11 piezoelectric patches equally spaced over a portion of the length. The wave propagation performance of the beam assembly is evaluated to illustrate the occurrence of attenuation bands over selected frequency ranges. The center frequencies for the attenuation bands are selected by tuning the electric impedance of the shunting circuit, which does not require any structural modification. Therefore, the proposed concept is an example of a waveguide with tunable attenuation properties.

Fig. 36
Fig. 36
Close modal

The resonant shunted circuits are implemented through the application of a synthetic inductor (Antoniou's circuit) because of the high value of inductance needed for the desired frequency tunings and for the flexibility in the circuit tuning procedure [399]. In the configuration tested, each pair of colocated piezoelectric patches is connected in series to a resonant shunting circuit. The electric networks were tuned at the frequencies of 5000 and 11,000 Hz. The beam is excited by the piezo actuator shown in Fig. 36, which is fed a pulse signal for broadband excitation of the beam motion. The velocity of the beam is measured over the entire beam length by the SLDV. Figure 37(a) presents an example of a recorded space-time plot, which shows two waves emanating from the excitation location at approximately x = 1 m, and subsequently reaching the ends of the beam where they are reflected.

Fig. 37
Fig. 37
Close modal

The estimation of the dispersion properties for the beam is performed by transforming the recorded response w(x, t) into the frequency/wavenumber through a 2D FT. This operation requires the preliminary removal of the boundary reflections in order to obtain clear dispersion representations. This is achieved by applying a dispersion compensation procedure that allows the identification of reflections as compact echoes and their elimination through windowing. The dispersion compensation is based on compensation conducted by applying the warped frequency transform (WFT). The result of this procedure leads to the space-time signal shown in Fig. 37(b), where the left boundary reflection is effectively removed. The windowed response is finally analyzed through a 2D FT, whose amplitude $|W(k,ω)|$ is represented as a surface in the frequency/wavenumber domain. The contour plot of the 2D FT surface for a beam with open circuits is shown in Fig. 38(a), where the maxima of the contours outline the dispersion relation for the beam. Within the considered frequency range, a continuous frequency/wavenumber relation can be observed for the beam with open circuits. In contrast, results for shunts tuned at 5000 Hz presented in Fig. 39 clearly show the presence of a gap in the frequency/wavenumber spectrum at the tuning frequency.

Fig. 38
Fig. 38
Close modal
Fig. 39
Fig. 39
Close modal
The spatial/temporal resolution provided by the considered experimental setup allows the quantitative estimation of wavenumbers from experimental data. Both real and imaginary part of the wavenumber are estimated as a function of frequency, so that attenuation frequency bands can be quantified. The approach is based on the estimation of the spatial variation of the FT of the response evaluated at each frequency of interest. The procedure is illustrated by considering the measured response as the superposition of harmonic waves of the kind $w(x,t)≈∑w∧i(x,t)$, where $w∧i(x,t)i=w∧0(ωi)ei[κ(ωi)x-ωit]$ is the ith harmonic component. The FT of $w∧i(x,t)i$ at ωi is approximately given by
$Wi(x,ωi)≈w∧0(ωi)eiκix$
(83)
where κi = κ(ωi). The wavenumber $κi=κiR+iκiI$ is generally a complex number, so that Eq. (83) can be rewritten as
$Wi(x,ωi)≈w∧0(ωi)e-κiIxeiκiRx$
(84)
Analysis of Eq. (84) reveals that the estimation of the FT of the recorded data and the evaluation of the variation of the resulting expression along the spatial coordinate x allow for the estimation of the spatial decay of the response amplitude and the phase evolution in space at a given frequency. This in turn leads to estimation of the attenuation constant as defined by the imaginary component of the wavenumber $κiI$, and of its real component, which are, respectively, given by
$κiI(xf-xi)=log[|W(x,ωi)|]$
(85)
$κiR(xf-xi)=arg[|W(x,ωi)|]$
(86)

where $x∈[xi xf]$, with xi, xf respectively denoting the initial coordinate and final coordinate over which the spatial decay and phase linear modulation are interpolated. Application of Eqs. (85) and (86) to the experimental results allows the quantification of the wavenumber variation over the frequency range of interest, and particularly the evaluation of the attenuation constant, examples of which are shown in Fig. 39. The red dashed lines correspond to the case of open circuits, while the solid blue lines are results for shunting at the tuning frequencies considered. These results clearly illustrate how shunting of the circuits at a given frequency creates an attenuation frequency band centered at the tuning frequency. Such a band is defined by large nonzero values of the imaginary part of the wavenumber, which is a measure of attenuation. In the same frequency range, the real part of the wavenumber defining the propagation component undergoes a resonant behavior consistent with theoretical results for internally resonating metamaterials. Of note is the fact that the case of open circuits does not lead to an absolute zero for the attenuation constant, which may be affected by other sources of dissipation that are inevitably present in an experimental setup.

### Stubbed Plate With Periodic Electromechanical Resonators.

The last example combines the two concepts presented above, i.e., a stubbed plate with an array of electromechanical resonators. The stubs in this case are all metallic and are designed to produce a Bragg scattering bandgap [213,214,400–402]. A 1D waveguide is designed in the stub array to localize wave motion within the bandgap frequency range. The addition of the piezoelectric resonators within the 1D waveguide generates a tunable resonance gap that can stop propagation of waves at the tuning frequencies. This example therefore documents the occurrence of a Bragg bandgap and the tunability of the waveguide propagation properties through internal resonance [258].

The surface-bonded cylindrical stubs have diameter d = 7 mm and height h = 10 mm. The stubs are made of aluminum and are arranged in a square lattice of constant a = 10 mm. An L-shaped waveguide is realized within a 17 × 17 lattice by removing parts of one row and one column of stubs as shown in Fig. 40. The vertical portion of the channel is equipped with a periodic array of 8 × 1 piezoceramic disks independently shunted through an inductive circuit.

Fig. 40
Fig. 40
Close modal

Wavefield measurements are conducted over a 2D grid by the SLDV. The waves are excited by a piezoelectric transducer located near the bend in the L-shaped waveguide (depicted in red in Fig. 40). As in the previous example, the measured data are analyzed through spatial and temporal FTs to eliminate disturbances from boundary reflections and to obtain direct estimates of real and imaginary parts of the wavenumber. The considered plate, originally proposed in Ref. [214], is characterized by a Bragg-type bandgap between 100 and 130 kHz [403]. This frequency range is defined by the material and geometry of plate and stubs and cannot be changed after the system is assembled. Figure 41(a) shows that the system behaves as an effective waveguide at frequencies inside the bandgap region illustrated in Fig. 41(b).

Fig. 41
Fig. 41
Close modal

The 2D map of the wavefield recorded with open circuits shown in Fig. 41(a) clearly illustrates how the wave transmission within the channel is not altered by the presence of the piezo disks. Also, Fig. 42 shows the contour of the frequency/wavenumber domain representation of the wavefield recorded along the vertical waveguide. The contours are overlapped to the theoretical dispersion branch of the first antisymmetric (A0) Lamb wave mode for the considered plate. The close match between the contours and the theoretical dispersion, and the continuous frequency/wavenumber distribution defined by the contours, suggest how the presence of the piezos with open circuits has little effect on wave transmission.

Fig. 42
Fig. 42
Close modal

Shunting of the piezos through the electrical networks produces a resonant, complex modulus (Eq. (82)), which in turn generates a frequency region of strong attenuation centered at ωT. The tunable properties of the waveguide are illustrated by wavefield measurements with circuits designed to resonate at different frequencies. Specifically, two sets of inductors (L1 = 1890μ H, L2 = 1560μ H) are used to tune the resonating units within the bandgap of the stubbed plate. In each case, contours of the frequency/wavenumber representation of the response, shown in Figs. 43(a) and 43(d), respectively, indicate that the resulting medium features a bandgap centered at the tuning frequencies. The measured data also allow the experimental evaluation of the attenuation and propagation constants. The real component of the wavenumber, overlapped as a solid line with square markers to the contour plots in Figs. 43(a) and 43(d), reveals the back-bending of the dispersion curve typical of internally resonating metamaterials. Also, the imaginary part, shown in Figs. 43(b) and 43(e), exhibits a sharp peak of attenuation in a small range of frequencies centered at the resonance of the electrical circuit ωT. The resulting attenuation properties of the waveguide are responsible for the onset of the tunable bandgap, and can be correlated to the frequency variation of the effective modulus of the resonating units as shown in Figs. 43(c) and 43(e).

Fig. 43
Fig. 43
Close modal

## Outlook

As documented herein, periodic structures, phononic crystals, and acoustic/elastic metamaterials are providing promising avenues for new and ongoing research, as well as a strong potential for transition of many concepts and fundamental studies to fully developed industrial technologies. The field is experiencing continued growth in research activities as demonstrated by the increasing number of symposia, workshops, and conferences,3 as well as journal special issues [286,404–407], devoted completely or in part to phononic materials and structures.

There are numerous directions in the study of dynamics of phononic materials and structures that we anticipate will occupy the imagination of researchers in the years to come. Seemingly unlimited future work promises to arise as a fruition of basic phenomena and analysis techniques developed in the past one to two decades (of which a substantial portion has been reviewed in this paper). Furthermore, it may still be possible to continue to discover phenomena in classical waves that bear resemblance to quantum waves, as was done with the observation of tunneling in phononic crystals, for example, Ref. [408].

Looking forward, we take, as a starting point, the fundamental topic of incorporation of damping and nonlinearity in phononic system models (Secs. 4.1 and 4.2, respectively). More complexity in the type of models, for example, an N-network generalized Maxwell viscoelastic model with nonlinear components, is yet to be developed. And equally as important is the investigation of experimental approaches for validating recent advancements in understanding dissipative and nonlinear Bloch wave propagation. While experiments have been conducted to probe the effects of damping and nonlinearity in finite structures, the extension to wave motion in periodic media is still at its infancy.

Another branch of research, also reviewed earlier in this article, is the study of coupling of elastic waves with other types of waves (especially electromagnetic). Much effort has been devoted to the exploration of unit-cell configurations that exhibit simultaneous phononic and photonic band gaps, and several experimental studies have successfully realized nonlinear multifield wave interactions. However, there is still room for further development of theoretical frameworks for modeling the coupled fields and their simultaneous wave propagation characteristics. This difficulty stems primarily from nonlinearities that may be associated with the coupling. Therefore, further advances in self-consistent nonlinear wave propagation analysis for periodic materials and structures will undoubtedly support this topic. With such development, the topics of phoXonic crystals and optomechanics will take strides toward the realization of devices that could serve numerous applications in both the fields of acoustics and optics.

Tunability is yet another area expected to continue to attract interest in the years to come as it represents a practical prerequisite to effective use of phononic materials and structures in engineering applications. While to date, the focus has been mainly on demonstration of different approaches and concepts for active control (as briefly surveyed in Sec. 1.3.2), future research will address the need to quantify the accuracy and speed at which band structure tuning may be conducted, and the subsequent rate at which the desired transformation of wave motion takes place. To this end, ideas and techniques from the well-established discipline of automatic control will be sought and applied to phononic systems.

The above categories of future research will collectively provide new opportunities for current and yet-to-be-conceived applications. For example, the promise of realization of phononic logic circuits and energy harvesting devices will become more realistic and will in turn bring about new opportunities. New applications and topics, such as dual control of vibrations and sound, incorporation of fluid-structure interactions and use of elastic phonons to control the buckling resistance of structured materials, may be also pursued.

Ultimately, the field will approach its two extreme frontiers: phononic control at the very low-frequency regime and at the very high-frequency regime. While acoustic/elastic metamaterials have provided a route to subwavelength applications, it is still a challenge to take the field to the ultralow-frequency domain. Should this be achieved, phononic materials and structures will impact a whole new category of applications ranging from vibration suppression of rotating machinery to earthquake mitigation, among others. The key challenge here is to realize subwavelength unit cells that are still capable of providing the desired dynamic performance, while at the same time exhibiting practical qualities, such as being small, compact, lightweight, and of practical implementation. At the ultrahigh-frequency extreme, the challenge is to take the current work on nanoscale phononic crystals to the point of full elucidation of the underlying interplay between coherent and incoherent scattering mechanisms. For this to be achieved, a key obstacle is the limited computational resources for adequate treatment of atomic wave motion at length scales spanning a few hundred nanometers. This suggests that advancements in unit-cell model reduction techniques at the atomic level are highly needed, and possible other ideas from multiscale modeling and coupling of continuum and atomistic descriptions. With the availability of these computational tools, as well as corresponding nanostructure fabrication and characterization technology, further crossfertilization will take place between concepts at the macro and nanoscales, such as the recently proposed notion of a “locally resonant nanophononic metamaterial.”

Should the current growth trend in phononic materials and structures research continue, it is not inconceivable for the field to witness a near-explosion in industrial interest. Such an outcome was witnessed over half a century ago in the field of electronics and more recently in the closer field of photonics. Phononics appear to be next.

## Acknowledgment

M. I. Hussein would like to acknowledge support by DARPA through its Young Faculty Award Grant YFA N66001-11-1-4134, support from the National Science Foundation through its Grant Nos. 0927322 and 1131802 and CAREER Grant No. 1254931, and support by the Defense Acquisition Program Administration and Agency for Defense Development under the Contract No. UD110096DD (Republic of Korea). He is also grateful to R. Khajehtourian, M. J. Frazier, and the rest of the CU-Boulder group for their assistance in proofreading the manuscript. M. J. Leamy and M. Ruzzene would like to acknowledge support from the National Science Foundation under Grant No. 0926776.

2

To the best of the authors' knowledge, this term was first used in the title of a symposium organized by J. S. Jensen at the 2005 U.S. National Congress in Computational Mechanics (Austin, TX, July 24–28, 2005).

3

Examples include the Annual ASME IMECE Symposium on Phononic Crystals and Acoustic Metamaterials (which started in 2005; http://www.asmeconferences.org/congress20xx), the International Workshop on Phononic Crystals (Nice, France, June 24–26, 2009; iwpc.gatech.edu) and the biannual Phononics 20xx conference series (whose inaugural edition was Phononics 2011: The First International Conference on Phononic Crystals, Metamaterials and Optomechanics, Santa Fe, NM, May 29–June 2, 2011; www.phononics2011.org).

## References

1.
Newton
,
I.
,
1686
,
Principia—Book II
,
Imprimatur
S.
Pepys
,
Reg. Soc. Præses
, London.
2.
Rayleigh
,
J. W. S.
,
1945
,
The Theory of Sound
, Vol.
1
,
Dover
,
New York
.
3.
Sun
,
C. T.
,
Achenbach
,
J. D.
, and
Herrmann
,
G.
,
1968
, “
Time-Harmonic Waves in a Stratified Medium Propagating in the Direction of the Layering
,”
ASME J. Appl. Mech.
,
35
(2), pp. 408–411.10.1115/1.3601212
4.
Nemat-Nasser
,
S.
,
1972
, “
General Variational Methods for Waves in Elastic Composites
,”
J. Elasticity
,
2
(
2
), pp.
73
90
.10.1007/BF00046056
5.
Abrahamson
,
A. L.
,
1973
, “
The Response of Periodic Structures to Aero-Acoustic Pressures With Particular Reference to Aircraft Skin-Rib-Spar Structures
,” Ph.D. thesis, University of Southampton, Southampton, U.K.
6.
,
D. J.
,
1973
, “
A General Theory of Harmonic Wave Propagation in Linear Periodic Systems With Multiple Coupling
,”
J. Sound Vib.
,
27
(
2
), pp.
235
260
.10.1016/0022-460X(73)90064-3
7.
Ewins
,
D. J.
,
1973
, “
Vibration Characteristics of Bladed Disk Assemblies
,”
J. Mech. Eng. Sci.
,
15
, pp.
165
186
.10.1243/JMES_JOUR_1973_015_032_02
8.
Griffin
,
J. H.
, and
Hoosac
,
T. M.
,
1984
, “
Model Development and Statistical Investigation of Turbine Blade Mistuning
,”
J. Vib., Acoust., Stress Reliab. Des.
,
106
, pp.
204
210
.10.1115/1.3269170
9.
Castanier
,
M. P.
,
,
G.
, and
Pierre
,
C.
,
1997
, “
A Reduced Order Modeling Technique for Mistuned Bladed Disks
,”
ASME J. Vib. Acoust.
,
119
(3), pp.
439
447
.10.1115/1.2889743
10.
Deshpande
,
V. S.
, and
Fleck
,
N. A.
,
2000
, “
High Strain Rate Compressive Behaviour of Aluminium Alloy Foams
,”
Int. J. Impact Eng.
,
24
, pp.
277
298
.10.1016/S0734-743X(99)00153-0
11.
Hazizan
,
M. A.
, and
Cantwell
,
W. J.
,
2002
, “
The Low Velocity Impact Response of Foam-Based Sandwich Structures
,”
Composites, Part B
,
33
, pp.
193
204
.10.1016/S1359-8368(02)00009-4
12.
Fleck
,
N. A.
, and
Deshpande
,
V. S.
,
2004
, “
,”
ASME J. Appl. Mech.
,
71
(3), pp.
386
401
.10.1115/1.1629109
13.
Talbot
,
J. P.
, and
Hunt
,
H. E. M.
,
2003
, “
A Computationally Efficient Piled-Foundation Model for Studying the Effects of Ground-Borne Vibration on Buildings
,”
Proc. Inst. Mech. Eng., Part C
,
217
, pp.
975
989
.10.1243/095440603322407227
14.
Bao
,
J.
,
Shi
,
Z. F.
, and
Xiang
,
H. J.
,
2012
, “
Dynamic Responses of a Structure With Periodic Foundations
,”
J. Eng. Mech.
,
138
, pp.
761
769
.10.1061/(ASCE)EM.1943-7889.0000383
15.
Brun
,
M.
,
Movchan
,
A. B.
, and
Jones
,
I. S.
,
2013
, “
Phononic Band Gap Systems in Structural Mechanics: Finite Slender Elastic Structures and Infinite Periodic Waveguides
,”
ASME J. Vib. Acoust.
,
135
(4), p.
041013
.10.1115/1.4023819
16.
Sigalas
,
M. M.
, and
Economou
,
E. N.
,
1992
, “
Elastic and Acoustic Wave Band Structure
,”
J. Sound Vib.
,
158
(
2
), pp.
377
382
.10.1016/0022-460X(92)90059-7
17.
Sigalas
,
M.
, and
Economou
,
E. N.
,
1993
, “
Band Structure of Elastic Waves in Two Dimensional Systems
,”
Solid State Commun.
,
86
, pp.
141
143
.10.1016/0038-1098(93)90888-T
18.
Kushwaha
,
M. S.
,
Halevi
,
P.
,
Dobrzynski
,
L.
, and
Djafari-Rouhani
,
B.
,
1993
, “
Acoustic Band Structure of Periodic Elastic Composites
,”
Phys. Rev. Lett.
,
71
(
13
), pp.
2022
2025
.10.1103/PhysRevLett.71.2022
19.
Kushwaha
,
M. S.
,
Halevi
,
P.
,
Martínez
,
G.
,
Dobrzynski
,
L.
, and
Djafari-Rouhani
,
B.
,
1994
, “
Theory of Acoustic Band Structure of Periodic Elastic Composites
,”
Phys. Rev. B
,
49
(
4
), pp.
2313
2322
.10.1103/PhysRevB.49.2313
20.
Vasseur
,
J. O.
,
Djafari-Rouhani
,
B.
,
Dobrzynski
,
L.
,
Kushwaha
,
M. S.
, and
Halevi
,
P.
,
1994
, “
Complete Acoustic Band Gaps in Periodic Fibre Reinforced Composite Materials: The Carbodepoxy Composite and Some Metallic Systems
,”
J. Phys.: Condens. Matter
,
6
, pp.
8759
8770
.
21.
Kushwaha
,
M. S.
,
1996
, “
Classical Band Structure of Periodic Elastic Composites
,”
Int. J. Mod. Phys. B
,
10
, pp.
977
1094
.10.1142/S0217979296000398
22.
Liu
,
Z. Y.
,
Zhang
,
X. X.
,
Mao
,
Y. W.
,
Zhu
,
Y. Y.
,
Yang
,
Z. Y.
,
Chan
,
C. T.
, and
Sheng
,
P.
,
2000
, “
Locally Resonant Sonic Materials
,”
Science
,
289
(
5485
), pp.
1734
1736
.10.1126/science.289.5485.1734
23.
Kittel
,
C.
, and
Kittel
,
B. C.
,
1976
,
Introduction to Solid State Physics
,
Wiley
,
New York
.
24.
Brillouin
,
L.
,
1946
,
Wave Propagation in Periodic Structures Electric Filters and Crystal Lattices
,
Dover Publications, Inc.
,
New York
.
25.
Kelvin
,
W. T.
,
1881
,
Popular Lectures
, Vol.
1
,
Macmillan and Company
, London.
26.
Rytov
,
S. M.
,
1956
, “
Acoustical Properties of a Thinly Laminated Medium
,”
Sov. Phys.-Acoust.
,
2
, pp.
68
80
.
27.
Sun
,
C. T.
,
Achenbach
,
J. D.
, and
Herrmann
,
G.
,
1968
, “
Continuum Theory for a Laminated Medium
,”
ASME J. Appl. Mech.
,
35
(3), p. 467–475.10.1115/1.3601237
28.
Nemat-Nasser
,
S.
,
1972
, “
Harmonic Waves in Layered Composites
,”
ASME J. Appl. Mech.
,
39
(3), pp. 850–852.10.1115/1.3422814
29.
Nemat-Nasser
,
S.
, and
Fu
,
F. C. L.
,
1974
, “
Harmonic Waves in Layered Composites: Bounds on Frequencies
,”
ASME J. Appl. Mech.
,
41
(1), pp. 288–290.10.1115/1.3423245
30.
Nemat-Nasser
,
S.
, and
Minagawa
,
S.
,
1975
, “
Harmonic Waves in Layered Composites: Comparison Among Several Schemes
,”
ASME J. Appl. Mech.
,
42
(3), pp. 699–704.10.1115/1.3423665
31.
Nemat-Nasser
,
S.
,
Fu
,
F. C. L.
, and
Minagawa
,
S.
,
1975
, “
Harmonic Waves in One-, Two- and Three-Dimensional Composites: Bounds for Eigenfrequencies
,”
Int. J. Solids Struct.
,
11
(
5
), pp.
617
642
.10.1016/0020-7683(75)90034-7
32.
Nemat-Nasser
,
S.
, and
,
M.
,
1981
, “
Harmonic Waves in Layered Transversely Isotropic Composites
,”
J. Sound Vib.
,
79
(
2
), pp.
161
170
.10.1016/0022-460X(81)90365-5
33.
Lee
,
E. H.
,
1972
, “
A Survey of Variational Methods for Elastic Wave Propagation Analysis in Composites With Periodic Structures
,”
Dynamics of Composite Materials
,
E. H.
Lee
, ed.,
ASME
,
New York
, pp.
122
138
.
34.
Hegemier
,
G. A.
, and
Nayfeh
,
A. H.
,
1973
, “
A Continuum Theory for Wave Propagation in Laminated Composites—Case 1: Propagation Normal to the Laminates
,”
ASME J. Appl. Mech.
,
40
, p.
503
.10.1115/1.3423013
35.
Nayfeh
,
A. H.
,
1974
, “
Time-Harmonic Waves Propagating Normal to the Layers of Multilayered Periodic Media
,”
ASME J. Appl. Mech.
,
41
(1), pp. 92–96.10.1115/1.3423282
36.
Nayfeh
,
A. H.
,
1991
, “
The General Problem of Elastic Wave Propagation in Multilayered Anisotropic Media
,”
J. Acoust. Soc. Am.
,
89
(
4
), pp.
1521
1531
.10.1121/1.400988
37.
Chimenti
,
D. E.
, and
Martin
,
R. W.
,
1991
, “
Nondestructive Evaluation of Composite Laminates by Leaky Lamb Waves
,”
Ultrasonics
,
29
(
1
), pp.
13
21
.10.1016/0041-624X(91)90168-8
38.
Shull
,
P. J.
,
Chimenti
,
D. E.
, and
Datta
,
S. K.
,
1994
, “
Elastic Guided Waves and the Floquet Concept in Periodically Layered Plates
,”
J. Acoust. Soc. Am.
,
95
, p.
99
.10.1121/1.408270
39.
Safaeinili
,
A.
, and
Chimenti
,
D. E.
,
1995
, “
Floquet Analysis of Guided Waves in Periodically Layered Composites
,”
J. Acoust. Soc. Am.
,
98
, p.
2336
.10.1121/1.413281
40.
Murakami
,
H.
, and
Akiyama
,
A.
,
1985
, “
A Mixture Theory for Wave Propagation in Angle-Ply Laminates, Part 2: Application
,”
ASME J. Appl. Mech.
,
52
, p.
338
.10.1115/1.3169050
41.
Murakami
,
H.
,
1985
, “
A Mixture Theory for Wave Propagation in Angle-Ply Laminates, Part 1: Theory
,”
ASME J. Appl. Mech.
,
52
(2), pp. 338–344.10.1115/1.3169049
42.
McDevitt
,
T. W.
,
Hulbert
,
G. M.
, and
Kikuchi
,
N.
,
1999
, “
Plane Harmonic Wave Propagation in Three-Dimensional Composite Media
,”
Finite Elem. Analysis Des.
,
33
(
4
), pp.
263
282
.10.1016/S0168-874X(99)00049-9
43.
McDevitt
,
T. W.
,
Hulbert
,
G. M.
, and
Kikuchi
,
N.
,
2001
, “
An Assumed Strain Method for the Dispersive Global-Local Modeling of Periodic Structures
,”
Computer Methods Appl. Mech. Eng.
,
190
(
48
), pp.
6425
6440
.10.1016/S0045-7825(00)00184-5
44.
Hussein
,
M. I.
,
Hulbert
,
G. M.
, and
Scott
,
R. A.
,
2006
, “
Mode-Enriched Dispersion Models of Periodic Materials Within a Multiscale Mixed Finite Element Framework
,”
Finite Elem. Anal. Des.
,
42
, pp.
602
612
.10.1016/j.finel.2005.11.002
45.
Fish
,
J.
, and
Belsky
,
V.
,
1995
, “
Multi-Grid Method for Periodic Heterogeneous Media Part 2: Multiscale Modeling and Quality Control in Multidimensional Case
,”
Computer Methods Appl. Mech. Eng.
,
126
(
1–2
), pp.
17
38
.10.1016/0045-7825(95)00812-F
46.
Fish
,
J.
,
Chen
,
W.
, and
Nagai
,
G.
,
2002
, “
Non-Local Dispersive Model for Wave Propagation in Heterogeneous Media: Multi-Dimensional Case
,”
Int. J. Numer. Methods Eng.
,
54
(
3
), pp.
347
363
.10.1002/nme.424
47.
Fish
,
J.
,
Chen
,
W.
, and
Nagai
,
G.
,
2002
, “
Non-Local Dispersive Model for Wave Propagation in Heterogeneous Media: One-Dimensional Case
,”
Int. J. Numer. Methods Eng.
,
54
(
3
), pp.
331
346
.10.1002/nme.423
48.
Cremer
,
L.
, and
Leilich
,
H. O.
,
1953
, “
Zur theorie der biegekettenleiter
,”
(On Theory of Flexural Periodic Systems), Arch. Elektr. Uebertrag
,
7
, pp.
261
270
.
49.
Müller
,
H. L.
,
1957
, “
Biegewellen-Dämmung an Symmetrischen und Exzentrischen Sperrmassen
,”
Frequenz
,
11
(
10
), pp.
325
331
.
50.
Heckl
,
M.
,
1961
, “
Wave Propagation on Beam-Plate Systems
,”
J. Acoust. Soc. Am.
,
33
, pp. 640–651.10.1121/1.1908750
51.
Heckl
,
M. A.
,
1964
, “
Investigations on the Vibrations of Grillages and Other Simple Beam Structures
,”
J. Acoust. Soc. Am.
,
36
, pp. 1335–1343.10.1121/1.1919206
52.
Ungar
,
E. E.
,
1966
, “
Steady-State Responses of One-Dimensional Periodic Flexural Systems
,”
J. Acoust. Soc. Am.
,
39
, pp. 887–894.10.1121/1.1909967
53.
Lin
,
Y. K.
, and
McDaniel
,
T. J.
,
1969
, “
Dynamics of Beam-Type Periodic Structures
,”
J. Eng. Ind.
,
91
, pp. 1133–1141.10.1115/1.3591761
54.
Faulkner
,
M. G.
, and
Hong
,
D. P.
,
1985
, “
Free Vibrations of a Mono-Coupled Periodic System
,”
J. Sound Vib.
,
99
, pp. 29–42.10.1016/0022-460X(85)90443-2
55.
McDaniel
,
T. J.
,
1971
, “
Dynamics of Circular Periodic Structures (Periodically Supported and Damped Closed Circular Beam Structure, Determining Frequency Response Matrix)
,”
J. Aircraft
,
8
, pp.
143
149
.10.2514/3.44245
56.
McDaniel
,
T. J.
, and
Chang
,
K. J.
,
1980
, “
Dynamics of Rotationally Periodic Large Space Structures
,”
J. Sound Vib.
,
68
(
3
), pp.
351
368
.10.1016/0022-460X(80)90392-2
57.
Thomas
,
D. L.
,
1979
, “
Dynamics of Rotationally Periodic Structures
,”
Int. J. Numer. Methods Eng.
,
14
(
1
), pp.
81
102
.10.1002/nme.1620140107
58.
Williams
,
F. W.
,
1986
, “
An Algorithm for Exact Eigenvalue Calculations for Rotationally Periodic Structures
,”
Int. J. Numer. Methods Eng.
,
23
(
4
), pp.
609
622
.10.1002/nme.1620230407
59.
Pierre
,
C.
,
1988
, “
Mode Localization and Eigenvalue Loci Veering Phenomena in Disordered Structures
,”
J. Sound Vib.
,
126
(
3
), pp.
485
502
.10.1016/0022-460X(88)90226-X
60.
Kim
,
M.
,
Moon
,
J.
, and
Wickert
,
J. A.
,
2000
, “
Spatial Modulation of Repeated Vibration Modes in Rotationally Periodic Structures
,”
ASME J. Vib. Acoust.
,
122
(1), pp. 62–68.10.1115/1.568443
61.
,
D. M.
,
1996
, “
Wave Propagation in Continuous Periodic Structures: Research Contributions From Southampton, 1964–1995
,”
J. Sound Vib.
,
190
(
3
), pp.
495
524
.10.1006/jsvi.1996.0076
62.
Gupta
,
S.
,
1970
, “
Dynamics of Periodically Stiffened Structures Using a Wave Approach
,” USAF Report No. AFML-TR-70-13.
63.
Gupta
,
S.
,
1970
, “
Natural Flexural Waves and the Normal Modes of Periodically-Supported Beams and Plates
,”
J. Sound Vib.
,
13
, pp.
89
101
.10.1016/S0022-460X(70)80082-7
64.
,
D. J.
,
1975
, “
Wave Propagation and Natural Modes in Periodic Systems: I. Mono-Coupled Systems
,”
J. Sound Vib.
,
40
, pp. 1–18.10.1016/S0022-460X(75)80227-6
65.
,
D. J.
,
1975
, “
Wave Propagation and Natural Modes in Periodic Systems: II. Multi-Coupled Systems, With and Without Damping
,”
J. Sound Vib.
,
40
(
1
), pp.
19
39
.10.1016/S0022-460X(75)80228-8
66.
,
D. J.
,
1986
, “
A New Method of Analyzing Wave Propagation in Periodic Structures; Applications to Periodic Timoshenko Beams and Stiffened Plates
,”
J. Sound Vib.
,
104
, pp. 9–27.10.1016/S0022-460X(86)80128-6
67.
,
D. J.
,
1970
, “
Vibration Response and Wave Propagation in Periodic Structures
,”
J. Manuf. Sci. Eng.
,
93
, pp. 783–792.10.1115/1.3428014
68.
,
D. J.
, and
Pujara
,
K. K.
,
1971
, “
Space Harmonic Analysis of Periodically Supported Beams: Response to Convected Random Loading
,”
J. Sound Vib.
,
14
, pp. 525–532.10.1016/0022-460X(71)90579-7
69.
,
D. J.
, and
Mallik
,
A. K.
,
1973
, “
An Approximate Method of Predicting the Response of Periodically Supported Beams Subjected to Random Convected Loading
,”
J. Sound Vib.
,
47
, pp. 457–471.10.1016/0022-460X(76)90873-7
70.
Orris
,
R. M.
, and
Petyt
,
M.
,
1974
, “
A Finite Element Study of Harmonic Wave Propagation in Periodic Structures
,”
J. Sound Vib.
,
33
(
2
), pp.
223
236
.10.1016/S0022-460X(74)80108-2
71.
Åberg
,
M.
, and
Gudmundson
,
P.
,
1997
, “
The Usage of Standard Finite Element Codes for Computation of Dispersion Relations in Materials With Periodic Microstructure
,”
J. Acoust. Soc. Am.
,
102
, pp. 2007–2013.10.1121/1.419652
72.
Mace
,
B. R.
,
Duhamel
,
D.
,
Brennan
,
M. J.
, and
Hinke
,
L.
,
2005
, “
Finite Element Prediction of Wave Motion in Structural Waveguides
,”
J. Acoust. Soc. Am.
,
117
, pp. 2835–2843.10.1121/1.1887126
73.
Manconi
,
E.
, and
Mace
,
B. R.
,
2009
, “
Wave Characterization of Cylindrical and Curved Panels Using a Finite Element Method
,”
J. Acoust. Soc. Am.
,
125
, pp. 154–163.10.1121/1.3021418
74.
Mace
,
B. R.
, and
Manconi
,
E.
,
2008
, “
Modelling Wave Propagation in Two-Dimensional Structures Using Finite Element Analysis
,”
J. Sound Vib.
,
318
(
4–5
), pp.
884
902
.10.1016/j.jsv.2008.04.039
75.
Langley
,
R. S.
,
1994
, “
On the Modal Density and Energy Flow Characteristics of Periodic Structures
,”
J. Sound Vib.
,
4
(
172
), pp.
491
511
.10.1006/jsvi.1994.1191
76.
Langley
,
R. S.
,
1996
, “
The Response of Two-Dimensional Periodic Structures to Point Harmonic Forcing
,”
J. Sound Vib.
,
197
, pp. 447–469.10.1006/jsvi.1996.0542
77.
Langley
,
R. S.
,
Bardell
,
N. S.
, and
Ruivo
,
H. M.
,
1997
, “
The Response of Two-Dimensional Periodic Structures to Harmonic Point Loading: A Theoretical and Experimental Study of Beam Grillage
,”
J. Sound Vib.
,
207
, pp. 521–535.10.1006/jsvi.1997.1154
78.
Ruzzene
,
M.
,
Mazzarella
,
L.
,
Tsopelas
,
P.
, and
Scarpa
,
F.
,
2002
, “
Wave Propagation in Sandwich Plates With Periodic Auxetic Core
,”
J. Intell. Mater. Syst. Struct.
,
13
(
9
), pp.
587
597
.10.1106/104538902031865
79.
Martinsson
,
P. G.
, and
Movchan
,
A. B.
,
2003
, “
Vibrations of Lattice Structures and Phononic Band Gaps
,”
Q. J. Mech. Appl. Math.
,
56
, pp.
45
64
.10.1093/qjmam/56.1.45
80.
Ruzzene
,
M.
,
Scarpa
,
F.
, and
Soranna
,
F.
,
2003
, “
Wave Beaming Effects in Two-Dimensional Cellular Structures
,”
Smart Mater. Struct.
,
12
, pp. 363–372.10.1088/0964-1726/12/3/307
81.
Ruzzene
,
M.
, and
Scarpa
,
F.
,
2005
, “
Directional and Band-Gap Behavior of Periodic Auxetic Lattices
,”
Phys. Status Solidi B
,
242
(
3
), pp.
665
680
.10.1002/pssb.200460385
82.
Phani
,
A. S.
,
Woodhouse
,
J.
, and
Fleck
,
N. A.
,
2006
, “
Wave Propagation in Two-Dimensional Periodic Lattices
,”
J. Acoust. Soc. Am.
,
119
, pp. 1995–2005.10.1121/1.2179748
83.
Gonella
,
S.
, and
Ruzzene
,
M.
,
2008
, “
Homogenization and Equivalent In-Plane Properties of Two-Dimensional Periodic Lattices
,”
Int. J. Solids Struct.
,
45
(
10
), pp.
2897
2915
.10.1016/j.ijsolstr.2008.01.002
84.
Wolfe
,
J.
,
1998
,
Imaging Phonons: Acoustic Wave Propagation in Solid.
,
Cambridge University Press
, Cambridge, U.K.
85.
Elachi
,
C.
,
1976
, “
Waves in Active and Passive Periodic Structures: A Review
,”
Proc. IEEE
,
64
, pp.
1666
1698
.10.1109/PROC.1976.10409