Non-reciprocal wave propagation in mechanically-modulated continuous elastic metamaterials

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,¹¹ 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,²⁶ and phononic crystals containing electromagnets.²⁷ 

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 testing^28,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,³⁵ 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.³⁶ 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) structure^31,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.³¹ 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.³⁸ 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.

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

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

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

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 as³⁹ 

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

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

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

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

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

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

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

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,

and for the linear elastic wave propagation

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

where β₀ is the static pre-strain, k_m = 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 L_X. 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.01 and Δβ = 0.3β₀ 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

where U = β₀L_y 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

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,

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

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

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.²⁰ 

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

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

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

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.⁴⁰ The supercell in Fig. 2(a) is spatially discretized into a computational mesh with a maximum element size of t_b/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.³¹ 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 L_ijkl 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

which is discretized following standard finite element procedures⁴¹ to yield a quadratic eigenvalue problem

where $U=[u−P,u−P+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, k_x is restricted to the first Brillouin zone of the supercell k_x ∈ [–π/(6L_x), π/(6 L_x)]. Positive values of k_x refer to waves that propagate in the positive x direction (with the modulation), and negative values of k_x 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 DOF_SC(2P + 1) × DOF_SC(2P + 1), where DOF_SC 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 SLEPc⁴² (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.¹⁸ The geometry is shown in Fig. 3, where h/L = 0.01. The plate is assigned artificial material properties of static Young's modulus E₀ = 1 GPa, density ρ₀ = 7700 kg/m³, 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

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

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

The modulation parameters are chosen to be E_m = 0.4E₀ and ω_m/k_m = 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.¹⁸ 

We note that this system is non-reciprocal since the frequency-wavenumber spectrum is not symmetric about k_x = 0. However, band gaps do not occur in the same way as a stationary periodic medium, where k_x 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 k_x 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 $20 log10(||û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/c₀ ∈ [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 ρ₀ = 1040 kg/m³, 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.³¹ with a static pre-strain of β₀ = 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.

Next, we obtain the frequency-wavenumber spectrum for modulation of the external pre-strain in space only by setting ω_m = 0 and Δβ = 0.3β₀ 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 c_m = 0.02c_s. 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).

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.¹⁸ 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.