Acoustic and elastic metamaterials with time- and space-dependent effective material properties have recently received significant attention as a means to induce non-reciprocal wave propagation. Recent analytical models of spring-mass chains have shown that external application of a nonlinear mechanical deformation, when applied on time scales that are slow compared to the characteristic times of propagating linear elastic waves, may induce non-reciprocity via changes in the apparent elastic modulus for perturbations around that deformation. Unfortunately, it is rarely possible to derive analogous analytical models for continuous elastic metamaterials due to complex unit cell geometry. The present work derives and implements a finite element approach to simulate elastic wave propagation in a mechanically-modulated metamaterial. This approach is implemented on a metamaterial supercell to account for the modulation wavelength. The small-on-large approximation is utilized to separate the nonlinear mechanical deformation (the “large” wave) from superimposed linear elastic waves (the “small” waves), which are then analyzed via Bloch wave analysis with a Fourier expansion in the harmonics of the modulation frequency. Results on non-reciprocal wave propagation in a negative stiffness chain, a structure exhibiting large stiffness modulations due to the presence of mechanical instabilities, are then shown as a case example.

For acoustic and elastic waves in linear-time-invariant heterogeneous media, the propagation of sound between two points is independent of the choice of source and receiver.1–3 This physical principle is known as acoustic reciprocity, and is generally obeyed except for certain specific scenarios. By breaking acoustic reciprocity, new avenues for direction-specific wave tailoring are opened, including the possibility of one-way sound propagation;2,4,5 this could assist in the design of direction-dependent acoustic devices with the potential to aid in numerous acoustical applications, such as energy harvesting, signal processing, vibration isolation, and acoustic communication.

Acoustic reciprocity can be broken by nonlinearity, which has been used to create one-way sound propagation via harmonic generation.6–10 Non-reciprocity can also be realized by applying a momentum bias, thereby breaking time-reversal symmetry, which has been achieved in piezophononic media,11 moving media,12,13 and gyroscopic phononic crystals,14,15 for example. A third mechanism, which is central to the present study, is modulation of the material properties in space and time.16–22 Effective mechanical properties have been manipulated in past works using electromagnetic effects, for example in piezoeletric materials,23–25 magnetorheological elastomers,26 and phononic crystals containing electromagnets.27 

One concept key to the present work is periodic, wave-like modulation in the form of nonlinear, purely mechanical deformation (the “large” wave), which effectively perturbs the linearized stiffness and/or mass properties of small disturbances propagating in superposition (the “small” wave). This “small-on-large” propagation behavior has been of interest for ultrasonic non-destructive testing28,29 and mechanical metamaterials.30–34 Small-on-large behavior has recently been used to achieve non-reciprocal wave propagation in chains of contacting cylinders,35 where the effective stiffness of longitudinal waves was modulated by an applied torsional deformation, as well as in another recent study, where the method of multiple scales was used to demonstrate spatiotemporal stiffness modulation in a periodic spring-mass lattice with a triangular unit cell.36 In both of these works, the considered structures could be modeled as quasi-one-dimensional lattices of lumped elements and/or rigid bodies, yielding analytical expressions for the modulated stiffness. However, in the case of more general structures, simple analytical models are not always available, and a more general formulation is needed.

In this work, we present a method for modeling continuous structures that achieve modulated effective stiffness via slow, nonlinear deformation. Using small-on-large theory, we formulate partial differential equations describing the nonlinear deformation and linear-elastic wave propagation. We then present a numerical implementation employing the finite element method, using a negative stiffness honeycomb (NSH) structure31,37 as a case example with complex geometry. The NSH is chosen because its incremental stiffness can be significantly altered by applying an external static pre-strain without operating near or past the buckling instabilities of the structure.31 The frequency-wavenumber spectrum of the non-reciprocal elastic waves is presented and compared with prior works. This method provides a framework for modeling non-reciprocal waves in spatiotemporally-modulated media when simple analytical models are not feasible.

In this section, we derive the small-on-large approximation for elastic waves in a general pre-stressed medium.38 The plane-strain approximation is assumed, reducing the material dynamics to two dimensions. In this work, we consider a quasi-one-dimensional chain of negative stiffness elements (termed negative stiffness chain hereafter), whose unit cell is depicted in Fig. 1. In order to capture the spatially varying pre-strain, a supercell must be modeled, which is constructed by repeating the unit cell in Fig. 1 along the x direction, as shown in Fig. 2(a). While the negative stiffness chain is considered here as a case example, this method can be applied to any material structure of interest.

FIG. 1.

Negative stiffness element with geometric parameters.

FIG. 1.

Negative stiffness element with geometric parameters.

Close modal
FIG. 2.

(Color online) (a) Supercell computational domain with displacement boundaries Γu and periodic boundary pairs (Γ+, Γ). (b) Nonlinear deformation of the negative stiffness chain due to the external pre-strain β(x, t) with parameters β0 = 0.01 and Δβ = 0.3β0. The external pre-strain is applied as displacement boundary conditions on the faces depicted by the arrows. The applied displacement in the y direction, uy(x, t), is composed of a static displacement term denoted by the dashed line whose value is β0Ly, and an oscillating component depicted by the sine wave with amplitude ΔβLy (not drawn to scale).

FIG. 2.

(Color online) (a) Supercell computational domain with displacement boundaries Γu and periodic boundary pairs (Γ+, Γ). (b) Nonlinear deformation of the negative stiffness chain due to the external pre-strain β(x, t) with parameters β0 = 0.01 and Δβ = 0.3β0. The external pre-strain is applied as displacement boundary conditions on the faces depicted by the arrows. The applied displacement in the y direction, uy(x, t), is composed of a static displacement term denoted by the dashed line whose value is β0Ly, and an oscillating component depicted by the sine wave with amplitude ΔβLy (not drawn to scale).

Close modal

Let ΩR2 be the material domain of a representative supercell in the undeformed configuration with coordinate X=[x,y]T. The equations of motion with respect to the reference configuration for the elastic medium are

DivS=ρ02ut2,
(1)

where u is the displacement, ρ0 is the density, and S is the first Piola-Kirchhoff stress tensor. S is related to the strain energy density function, W, by

S=WF,
(2)

where F is the deformation gradient, which is defined in relation to the displacement as

F=uX+I,
(3)

where I is the identity tensor. The particular strain energy density function is chosen to be the St. Venant-Kirchhoff model, which accounts for finite displacements but assumes small local strain values, and is written as39 

W=λ2tr(E)2+μtr(E2),
(4)
E=FTFI,
(5)

where λ and μ are the first and second Lamé parameters of the elastic medium, respectively, and E is the Green-Lagrangian strain tensor.

The total displacement u is assumed to be of the form

u(X,t)=ud(X,t)+ua(X,t),
(6)

where ud is the displacement due to the external pre-strain, and ua is the propagating elastic wave, which can be thought of as a perturbation, or incremental motion, about ud. The deformation gradient can thus be decomposed by substituting Eq. (6) into Eq. (3), which yields

F=Fd+Fa,
(7)

where Fd is the geometrically nonlinear deformation gradient due to the external pre-strain and is written as

Fd=udX+I,
(8)

and Fa is the deformation gradient due to the elastic wave motion, which is found to be

Fa=uaX.
(9)

The first Piola-Kirchhoff stress tensor can be approximated as a Taylor series expansion about the pre-strained state

SSd+Sa,
(10)

where Sd = S(Fd) is the stress due to the external pre-strain, and the incremental stress Sa can be written in index notation as

(Sa)ij=Lijkl(Fa)kl,
(11)

where Lijkl is the fourth-order tangent modulus tensor, which is defined as

Lijkl=2WFijFkl|u=ud.
(12)

Equation (1) can therefore be decomposed into the deformation and incremental terms using Eqs. (6) and (10) to yield a set of differential equations for the nonlinear deformation due to an external pre-strain,

DivSd=ρ02udt2,
(13)

and for the linear elastic wave propagation

Div(Lua)=ρ02uat2.
(14)

In the presence of spatiotemporally-modulated pre-strain, Eq. (14) is a wave equation with non-constant coefficients that are functions of space and time. The deformation due to an external pre-strain therefore acts as an effective stiffness change (in the reference configuration).

The displacement field due to an applied pre-strain is modeled using Eq. (13). The applied pre-strain, β, is assumed to have the following periodic traveling wave form

β(x,t)=Δβsin(kmxωmt)+β0,
(15)

where β0 is the static pre-strain, km = 2π/λm is the modulation wavenumber with wavelength λm, and ωm is the modulation angular frequency. For computational simplicity, we restrict λm to be equal to an integer number of unit cells. In this work, we consider a supercell composed of six unit cells, such that the supercell length and modulation wavelength are λm = 6 LX. The computational geometry is shown in Fig. 2(a), with boundaries Γu indicating where the external pre-strain is applied via displacement boundary conditions, and the boundary pairs (Γ+, Γ) indicating where periodic boundary conditions are applied. The displacements on the top and bottom faces are ud=(0,Lyβ) and ud=(0,Lyβ), respectively. The resulting deformation of the lattice due to the applied pre-strain with parameters β0 = 0.01 and Δβ = 0.3β0 at one instant in time is shown in Fig. 2(b). The depicted sine wave represents Eq. (15), and the arrows indicate where the displacements are applied. The external pre-strain then translates with time in the positive x direction.

The numerical modeling is greatly simplified and the displacement solutions are guaranteed to be periodic in time and space by requiring that the pre-strain modulation be slow with respect to the slowest sound speed in the medium (shear wave speed). To investigate the consequences of a slow modulation, it is illustrative to rewrite Eq. (13) as a non-dimensional equation. By choosing the normalization X¯=kmX,u¯d=ud/U,S¯d=Sd/(μUkm),t¯=ωmt, Eq. (13) is rewritten as

DivS¯d=γ2u¯dt¯2,
(16)

where U = β0Ly is a reference displacement, γ=cm2/cs2,cm=ωm/km is the modulation speed, and cs=μ/ρ0 is the shear wave speed of the material. If the modulation speed is much slower than the shear wave speed (γ ≪1), the inertial term can be discarded and Eq. (13) simplifies to

DivSd=0.
(17)

Solutions for ud(X,t) are constructed by discretizing t and incrementally solving Eq. (17) by updating the displacement boundary conditions. Consequently, nonlinear propagation effects such as harmonic generation or shock formation are absent.

Once the displacements ud are found, the tensor L is constructed for each material point and time increment using Eq. (12). Each element of L can then be written as a Fourier series in time,

Lijkl(X,t)=n=L̂ijkln(X)einωmt,
(18)

where the Fourier components, L̂ijkln(X), are determined from the set of solutions to Eq. (17) for all time increments of the modulation, using

L̂ijkln(X)=ωm2ππ/ωmπ/ωmLijkl(X,t)einωmtdt.
(19)

A Bloch wave solution is utilized to find the periodic traveling wave solutions of ua in Eq. (14) using the ansatz

(ua)i=ei(K·Xωt)p=ûip(X)eipωmt,
(20)

where K=[kx,ky]T is the Bloch wavenumber, and the Bloch wave mode is written as a Fourier series in time. Physically speaking, the fundamental mode and frequency (û0,ω) is interpreted as an incident mode propagating in the medium, and the harmonics with mode and frequency (ûp,ω+pωm) are modes scattered by the modulation in the forward (p >0) and backward (p <0) direction.20 

The series in Eqs. (18) and (20) are truncated to 2 P + 1 terms, and substituted into Eq. (14). The orthogonality of the Fourier series is then utilized to eliminate one of the summations, yielding a system of differential equations for the harmonic amplitudes ûp, which is written in index notation as

n=PP[(L̂ijklpnûk,ln+iL̂ijklpnKlûkn),j+iL̂ijklpnûk,lnKjL̂ijklpnûknKlKj]=ρ0(ω+pωm)2ûip,p[P,P].
(21)

This set of equations is then solved subject to homogeneous Dirichlet boundary conditions on the faces where the external pre-strain was applied,

ûp(XΓu)=0,
(22)

and periodic boundary conditions on the left and right edges as indicated in Fig. 2(a),

ûp(XΓ+)=ûp(XΓ).
(23)

The finite element method is used to solve both the nonlinear deformation and elastic wave equations derived in Sec. II by implementing the finite element method using the open source software FEniCS.40 The supercell in Fig. 2(a) is spatially discretized into a computational mesh with a maximum element size of tb/4 to accurately resolve the geometry of the negative stiffness chain and each Bloch wave mode at the frequencies of interest. The modulation deformation problem is solved for each time step using the same techniques outlined in Goldsberry et al.31 Once ud is computed, the tangent modulus tensor L is created using an automatic differentiation algorithm available in FEniCS. The Fourier components L̂ijkl are then found by taking the time solution of Lijkl at each node in the mesh, performing a numerical Fast Fourier Transform (FFT) in time, and reconstructing the results of the FFT into a set of tensors {L̂P,L̂P+1,,L̂P}.

In order to implement Eq. (21) in a finite element algorithm, the appropriate weak forms must be derived. This is accomplished by taking the Hermitian inner product over the supercell domain Ω of the left- and right-hand sides of Eq. (21) with the respective harmonic test vector v¯p. Application of Green's identity on the divergence term and summing the equations yields the integral equation

p,n=PP[ΩL̂ijklpn(ûk,lnv̂i,jp+KlKjûknv̂ip)dΩiΩL̂ijklpn(ûknKlv̂i,jpûk,lnKjv̂ip)dΩ+ρ0(ω+pωm)2Ωûipv̂ipdΩ]=0,
(24)

which is discretized following standard finite element procedures41 to yield a quadratic eigenvalue problem

ω2MU+ωCU+K(K)U=0,
(25)

where U=[uP,uP+1,,uP]T. The frequency-wavenumber spectrum is generated by assigning a value of K and solving the quadratic eigenvalue problem for the frequency and harmonic amplitudes. In this work, we restrict K=[kx,0]T to be a vector pointing in the x direction. Due to the periodicity of the supercell, kx is restricted to the first Brillouin zone of the supercell kx ∈ [–π/(6Lx), π/(6 Lx)]. Positive values of kx refer to waves that propagate in the positive x direction (with the modulation), and negative values of kx refer to waves that propagate in the negative x direction (against the modulation).

The matrices resulting from the finite element discretization are large and sparse. Specifically, the total matrix size is DOFSC(2P + 1) × DOFSC(2P + 1), where DOFSC is the number of degrees of freedom of the supercell with boundary conditions applied. Therefore, it is not practical to solve for all eigenfrequencies and eigenvectors. To address the computational challenges associated with this problem, we use the software library SLEPc42 (the Scalable Library for Eigenvalue Problem Computations) to solve Eq. (25), which takes advantage of distributed-memory parallelization and can be executed on a computer cluster. Further, a Two-Level Orthogonal Arnoldi (TOAR) algorithm is used in combination with a shift-and-invert transformation to extract eigenfrequencies near a target magnitude value.

The finite element approach described in Sec. III can be used to compute the frequency-wavenumber spectrum for any geometry of interest. Before solving the negative stiffness chain problem, we first assess the accuracy of the finite element model by comparing the frequency-wavenumber spectrum calculated with this technique to wave propagation in a thin Kirchhoff plate with a spatiotemporally-modulated Young's modulus, which has an analytical solution. We then demonstrate the use of the finite element model by investigating transverse wave propagation in the negative stiffness chain with a spatiotemporally-modulated external pre-strain.

To verify the implementation of the finite element model presented in Sec. III, we investigate elastic wave dispersion in a thin Kirchhoff plate with a spatiotemporally-modulated Young's modulus and compare with the model presented in Trainiti et al.18 The geometry is shown in Fig. 3, where h/L =0.01. The plate is assigned artificial material properties of static Young's modulus E0 = 1 GPa, density ρ0 = 7700 kg/m3, and Poisson's ratio ν = 0.1. These values are chosen such that Kirchhoff plate theory remains valid for all wavenumbers of interest. The Young's modulus is modulated with the form

E(x,t)=E0+Emcos(ωmtkmx).
(26)

For this case, the tangent modulus tensor in Eq. (12) simplifies to the stiffness tensor

Lijkl=λδijδkl+μ(δikδjl+δilδjk),
(27)

where δij is the Kronecker delta function, and λ, μ are the first and second Lamé parameters of the plate, respectively, which are related to the Young's modulus and Poisson's ratio through the relations

λ=Eν(1+ν)(12ν),μ=E2(1+ν).
(28)

The modulation parameters are chosen to be Em = 0.4E0 and ωm/km = 0.002. The transverse mode branch with P =1 is shown in Fig. 4(a), where c0=E0/ρ is the longitudinal wave speed of the plate. As seen in Fig. 4(a), the finite element model shows excellent agreement with the model presented in Trainiti et al.18 

FIG. 3.

Geometry of the thin plate Kirchhoff benchmark system.

FIG. 3.

Geometry of the thin plate Kirchhoff benchmark system.

Close modal
FIG. 4.

(a) (Color online) Frequency-wavenumber spectrum of the transverse wave in a thin Kirchhoff plate. Open circles are the results from Trainiti et al. (Ref. 18), and filled circles are results obtained from the finite element model. (b) The finite element results presented in (a) but with each point colored by the magnitude of the fundamental component in decibels, 20log10(||û0||/||U||).

FIG. 4.

(a) (Color online) Frequency-wavenumber spectrum of the transverse wave in a thin Kirchhoff plate. Open circles are the results from Trainiti et al. (Ref. 18), and filled circles are results obtained from the finite element model. (b) The finite element results presented in (a) but with each point colored by the magnitude of the fundamental component in decibels, 20log10(||û0||/||U||).

Close modal

We note that this system is non-reciprocal since the frequency-wavenumber spectrum is not symmetric about kx = 0. However, band gaps do not occur in the same way as a stationary periodic medium, where kx becomes purely imaginary in certain frequency regions. Rather, due to the existence of harmonics in a spatiotemporally-modulated medium, all frequencies are associated with at least one real-valued kx that represents a propagating mode. However, the amplitudes of the harmonic modes [p ≠ 0 in Eq. (20)] may be significantly less in magnitude than the fundamental mode (p =0 term). The amplitude differences are indicated in Fig. 4(b), where the finite element results in Fig. 4(a) are plotted and assigned a color that is determined by the magnitude of the fundamental component of the Bloch wave solution in decibels, specifically 20log10(||û0||/||U||). One property of the quadratic eigenvalue problem provided in Eq. (25) is that the eigenvectors for each distinct eigenvalue are no longer orthogonal. Physically, this means that each harmonic mode up can couple and exchange energy. Previous works have shown that directional band gaps occur when two or modes interact.18,20 In the example case considered here, the fundamental mode couples to the harmonic modes in the frequency range fL/c0 ∈ [0.003, 0.006]. However, this coupling is direction-dependent, creating directional band gaps. The amplitudes of the harmonic modes within the fundamental mode band gaps are significantly lower in magnitude than the amplitude of the fundamental mode propagating in the −x direction.

The negative stiffness chain depicted in Fig. 2(a) is now investigated in order to demonstrate the generality of the finite element method and to investigate mechanical modulation as a means to generate non-reciprocal wave phenomena. The values assigned to the geometric parameters shown in Fig. 1 are provided in Table I, and the material is made of laser sintered Nylon 11 with density ρ0 = 1040 kg/m3, Poisson's ratio ν = 0.33, and Young's modulus E =1582 MPa. The frequency-wavenumber spectrum of the negative stiffness element depicted in Fig. 1 is first studied using the finite element approach developed in Goldsberry et al.31 with a static pre-strain of β0 = 0.01. The frequency-wavenumber spectrum for the horizontal Brillouin zone is shown in Fig. 5(a). Only the transverse mode, which is highlighted in Fig. 5(a) and repeated with an enlarged view in Fig. 5(b), is investigated in the present work. The vertical dashed lines and numbers in Fig. 5(b) represent the Brillouin zones of the supercell shown in Fig. 2(a), which will aid in the interpretation of the band-folded supercell results detailed below.

TABLE I.

Geometric parameter values depicted in Fig. 1. All numerical values have units of mm.

ParameterDescriptionValue
Lx Horizontal length 55.88 
Ly Vertical length 40.64 
tb Beam thickness 1.27 
ts Beam separation 3.81 
hb Beam apex height 5.08 
hc Center height 1.90 
wc Center width 3.8 
hcb Center beam height 2.54 
wcb Center beam width 2.54 
thb Horizontal beam thickness 1.27 
ParameterDescriptionValue
Lx Horizontal length 55.88 
Ly Vertical length 40.64 
tb Beam thickness 1.27 
ts Beam separation 3.81 
hb Beam apex height 5.08 
hc Center height 1.90 
wc Center width 3.8 
hcb Center beam height 2.54 
wcb Center beam width 2.54 
thb Horizontal beam thickness 1.27 
FIG. 5.

(a) (Color online) Frequency-wavenumber spectrum of the negative stiffness chain with a static pre-strain of β0 = 0.01. The transverse mode is highlighted. (b) The transverse mode spectrum with the vertical lines and numerals indicating the supercell Brillouin zone number. The mode shape is shown on the top right.

FIG. 5.

(a) (Color online) Frequency-wavenumber spectrum of the negative stiffness chain with a static pre-strain of β0 = 0.01. The transverse mode is highlighted. (b) The transverse mode spectrum with the vertical lines and numerals indicating the supercell Brillouin zone number. The mode shape is shown on the top right.

Close modal

Next, we obtain the frequency-wavenumber spectrum for modulation of the external pre-strain in space only by setting ωm = 0 and Δβ = 0.3β0 in Eq. (15). The results of this case are shown in Fig. 6(a). The frequency-wavenumber spectrum of the transverse mode in a spatially-modulated negative stiffness chain is nearly identical to the folded frequency-wavenumber spectrum of the static pre-strain case shown in Fig. 5(b), with band gaps forming near 990 and 1410 Hz at the edge of the supercell Brillouin zone due to Bragg scattering from the spatial periodicity. Finally, the frequency-wavenumber spectrum of the transverse mode in the spatiotemporally-modulated pre-strain case with eight of the generated harmonics present (P =4) is shown in Fig. 6(b) with cm = 0.02cs. This modulation speed is chosen such that the small-on-large approximation is valid. The transverse mode fundamental branch is similar to the results from the spatially-modulated negative stiffness chain except at frequencies near the band gaps. Here, directional band gaps are now present in the spatiotemporally-modulated negative stiffness chain and are highlighted with shaded boxes in Fig. 6(b).

FIG. 6.

(Color online) Frequency-wavenumber spectra of the transverse mode of the negative stiffness chain. (a) Comparison of spatial modulation (open circles) and static pre-strain only [diamonds; same data as in Fig. 5(b), but folded at the Brillouin zone boundaries]. (b) Comparison of spatial modulation [open circles; same as in (a)] and spatiotemporal modulation (filled circles), where the color scale is defined in the same way as Fig. 4(b).

FIG. 6.

(Color online) Frequency-wavenumber spectra of the transverse mode of the negative stiffness chain. (a) Comparison of spatial modulation (open circles) and static pre-strain only [diamonds; same data as in Fig. 5(b), but folded at the Brillouin zone boundaries]. (b) Comparison of spatial modulation [open circles; same as in (a)] and spatiotemporal modulation (filled circles), where the color scale is defined in the same way as Fig. 4(b).

Close modal

We have developed a finite element model for non-reciprocal wave propagation in continuous elastic metamaterials, where the modulation is achieved by nonlinear deformation from an applied external pre-strain. We have assessed the accuracy of the finite element formulation by comparing the frequency-wavenumber spectrum of a thin Kirchhoff plate with an analytical model.18 Finally, we have demonstrated the procedure on a negative stiffness chain, which displays directional band gaps and mode coupling for the transverse wave modes.

This model can be used to study sub-wavelength geometries or modulations that are difficult or impossible to model with analytic techniques. Therefore, this approach can be used to design and optimize devices that benefit from a large degree of non-reciprocity, such as acoustic communication devices with increased data throughput and improved vibration isolation devices. Future work includes incorporating the finite element approach into a design optimization algorithm and investigating the degree of non-reciprocity as a function of the geometric parameters.

This work supported by National Science Foundation EFRI Award No. 1641078 and the postdoctoral fellowship program at Applied Research Laboratories at The University of Texas at Austin.

1.
J. W.
Strutt
, “
Some general theorems relating to vibrations
,”
P. Lond. Math Soc.
s1-4
(
1
),
357
368
(
1871
).
2.
R.
Fleury
,
D. L.
Sounas
,
C. F.
Sieck
,
M. R.
Haberman
, and
A.
Alù
, “
Sound isolation and giant linear nonreciprocity in a compact acoustic circulator
,”
Science
343
(
6170
),
516
519
(
2014
).
3.
J.
Achenbach
,
Reciprocity in Elastodynamics
(
Cambridge University Press
,
Cambridge, UK
,
2003
).
4.
M. R.
Haberman
and
M. D.
Guild
, “
Acoustic metamaterials
,”
Phys. Today
69
(
6
),
42
48
(
2016
).
5.
S. A.
Cummer
,
J.
Christensen
, and
A.
Alù
, “
Controlling sound with acoustic metamaterials
,”
Nat. Rev. Mater.
1
(
3
),
16001
(
2016
).
6.
N.
Boechler
,
G.
Theocharis
, and
C.
Daraio
, “
Bifurcation-based acoustic switching and rectification
,”
Nat. Mater.
10
(
9
),
665
668
(
2011
).
7.
Z.
Zhang
,
I.
Koroleva
,
L. I.
Manevitch
,
L. A.
Bergman
, and
A. F.
Vakakis
, “
Nonreciprocal acoustics and dynamics in the in-plane oscillations of a geometrically nonlinear lattice
,”
Phys. Rev. E
94
,
032214
(
2016
).
8.
J.
Bunyan
,
K. J.
Moore
,
A.
Mojahed
,
M. D.
Fronk
,
M.
Leamy
,
S.
Tawfick
, and
A. F.
Vakakis
, “
Acoustic nonreciprocity in a lattice incorporating nonlinearity, asymmetry, and internal scale hierarchy: Experimental study
,”
Phys. Rev. E
97
(
5
),
052211
(
2018
).
9.
B.
Liang
,
X.
Guo
,
J.
Tu
,
D.
Zhang
, and
J.
Cheng
, “
An acoustic rectifier
,”
Nat. Mater.
9
(
12
),
989
(
2010
).
10.
B.
Liang
,
B.
Yuan
, and
J.-C.
Cheng
, “
Acoustic diode: Rectification of acoustic energy flux in one-dimensional systems
,”
Phys. Rev. Lett.
103
(
10
),
104301
(
2009
).
11.
A.
Merkel
,
M.
Willatzen
, and
J.
Christensen
, “
Dynamic nonreciprocity in loss-compensated piezophononic media
,”
Phys. Rev. Appl.
9
(
3
),
034033
(
2018
).
12.
O. A.
Godin
, “
Reciprocity and energy theorems for waves in a compressible inhomogeneous moving fluid
,”
Wave Motion
25
(
2
),
143
167
(
1997
).
13.
O. A.
Godin
, “
Recovering the acoustic green's function from ambient noise cross correlation in an inhomogeneous moving medium
,”
Phys. Rev. Lett.
97
(
5
),
054301
(
2006
).
14.
L. M.
Nash
,
D.
Kleckner
,
A.
Read
,
V.
Vitelli
,
A. M.
Turner
, and
W. T. M.
Irvine
, “
Topological mechanics of gyroscopic metamaterials
,”
Proc. Natl. Acad. Sci. U.S.A.
112
(
47
),
14495
14500
(
2015
).
15.
P.
Wang
,
L.
Lu
, and
K.
Bertoldi
, “
Topological phononic crystals with one-way elastic edge waves
,”
Phys. Rev. Lett.
115
,
104302
(
2015
).
16.
E.
Cassedy
and
A.
Oliner
, “
Dispersion relations in time-space periodic media: Part I—Stable interactions
,”
Proc. IEEE
51
(
10
),
1342
1359
(
1963
).
17.
E.
Cassedy
, “
Dispersion relations in time-space periodic media Part II—Unstable interactions
,”
Proc. IEEE
55
(
7
),
1154
1168
(
1967
).
18.
G.
Trainiti
and
M.
Ruzzene
, “
Non-reciprocal elastic wave propagation in spatiotemporal periodic structures
,”
New J. Phys.
18
(
8
),
083047
(
2016
).
19.
J.
Vila
,
R. K.
Pal
,
M.
Ruzzene
, and
G.
Trainiti
, “
A bloch-based procedure for dispersion analysis of lattices with periodic time-varying properties
,”
J. Sound Vib.
406
,
363
377
(
2017
).
20.
H.
Nassar
,
H.
Chen
,
A. N.
Norris
,
M. R.
Haberman
, and
G. L.
Huang
, “
Non-reciprocal wave propagation in modulated elastic metamaterials
,”
Proc. Royal Soc. Lond. A
473
(
2202
),
20170188
(
2017
).
21.
H.
Nassar
,
X.
Xu
,
A. N.
Norris
, and
G. L.
Huang
, “
Modulated phononic crystals: Non-reciprocal wave propagation and Willis materials
,”
J. Mech. Phys. Solids
101
,
10
29
(
2017
).
22.
H.
Nassar
,
H.
Chen
,
A. N.
Norris
, and
G. L.
Huang
, “
Quantization of band tilting in modulated phononic crystals
,”
Phys. Rev. B
97
(
1
),
014305
(
2018
).
23.
F.
Casadei
,
T.
Delpero
,
A.
Bergamini
,
P.
Ermanni
, and
M.
Ruzzene
, “
Piezoelectric resonator arrays for tunable acoustic waveguides and metamaterials
,”
J. Appl. Phys.
112
(
6
),
064902
(
2012
).
24.
Y. Y.
Chen
,
G. L.
Huang
, and
C. T.
Sun
, “
Band gap control in an active elastic metamaterial with negative capacitance piezoelectric shunting
,”
J. Vib. Acoust.
136
(
6
),
061008
(
2014
).
25.
Y. Y.
Chen
,
R.
Zhu
,
M. V.
Barnhart
, and
G. L.
Huang
, “
Enhanced flexural wave sensing by adaptive gradient-index metamaterials
,”
Sci. Rep.
6
(
1
),
35048
(
2016
).
26.
K.
Danas
,
S.
Kankanala
, and
N.
Triantafyllidis
, “
Experiments and modeling of iron-particle-filled magnetorheological elastomers
,”
J. Mech. Phys. Solids
60
(
1
),
120
138
(
2012
).
27.
Y.
Wang
,
B.
Yousefzadeh
,
H.
Chen
,
H.
Nassar
,
G.
Huang
, and
C.
Daraio
, “
Observation of nonreciprocal wave propagation in a dynamic phononic lattice
,”
Phys. Rev. Lett.
121
,
194301
(
2018
).
28.
G.
Renaud
,
S.
Callé
, and
M.
Defontaine
, “
Remote dynamic acoustoelastic testing: Elastic and dissipative acoustic nonlinearities measured under hydrostatic tension and compression
,”
Appl. Phys. Lett.
94
(
1
),
011905
(
2009
).
29.
Y.
Zhang
,
V.
Tournat
,
O.
Abraham
,
O.
Durand
,
S.
Letourneur
,
A.
Le Duff
, and
B.
Lascoup
, “
Nonlinear mixing of ultrasonic coda waves with lower frequency-swept pump waves for a global detection of defects in multiple scattering media
,”
J. Appl. Phys.
113
(
6
),
064905
(
2013
).
30.
K.
Bertoldi
and
M. C.
Boyce
, “
Wave propagation and instabilities in monolithic and periodically structured elastomeric materials undergoing large deformations
,”
Phys. Rev. B
78
(
18
),
184107
(
2008
).
31.
B. M.
Goldsberry
and
M. R.
Haberman
, “
Negative stiffness honeycombs as tunable elastic metamaterials
,”
J. Appl. Phys.
123
(
9
),
091711
(
2018
).
32.
A.
Amendola
,
A.
Krushynska
,
C.
Daraio
,
N. M.
Pugno
, and
F.
Fraternali
, “
Tuning frequency band gaps of tensegrity mass-spring chains with local and global prestress
,”
Int. J. Solids Struct.
155
,
47
56
(
2018
).
33.
A. N.
Norris
and
W. J.
Parnell
, “
Hyperelastic cloaking theory: Transformation elasticity with pre-stressed solids
,”
Proc. R. Soc. A
468
(
2146
),
2881
2903
(
2012
).
34.
P.
Zhang
and
W. J.
Parnell
, “
Soft phononic crystals with deformation-independent band gaps
,”
Proc. R. Soc. A
473
(
2200
),
20160865
(
2017
).
35.
R.
Chaunsali
,
F.
Li
, and
J.
Yang
, “
Stress wave isolation by purely mechanical topological phononic crystals
,”
Sci. Rep.
6
(
1
),
30662
(
2016
).
36.
S. P.
Wallen
and
M. R.
Haberman
, “
Nonreciprocal wave phenomena in spring-mass chains with effective stiffness modulation induced by geometric nonlinearity
,”
Phys. Rev. E
99
,
013001
(
2019
).
37.
D. M.
Correa
,
T.
Klatt
,
S.
Cortes
,
M.
Haberman
,
D.
Kovar
, and
C.
Seepersad
, “
Negative stiffness honeycombs for recoverable shock isolation
,”
Rapid Prototyp. J.
21
(
2
),
193
200
(
2015
).
38.
R. W.
Ogden
, “
Incremental statics and dynamics of pre-stressed elastic materials
,” in
Waves in Nonlinear pre-Stressed Materials
(
Springer
,
New York
,
2007
), pp.
1
26
.
39.
K. D.
Hjelmstad
,
Fundamentals of Structural Mechanics
(
Springer Science & Business Media
,
New York
,
2007
).
40.
A.
Logg
,
K.-A.
Mardal
, and
G. N.
Wells
,
Automated Solution of Differential Equations by the Finite Element Method
(
Springer
,
New York
,
2012
).
41.
M. S.
Gockenbach
,
Understanding and Implementing the Finite Element Method
(
SIAM
,
Philadelphia
,
2006
).
42.
J. E.
Roman
,
C.
Campos
,
E.
Romero
, and
A.
Tomas
, “
SLEPc Users Manual
,” DSIC-II/24/02 - Revision 3.10, D (Sistemes Informàtics i Computació, Universitat Politècnica de València,
2018
).