Negative stiffness honeycombs as tunable elastic metamaterials

One of the principle objectives of elastic and acoustic metamaterial research is the design of heterogeneous material systems whose subwavelength structure generates macroscopically observable wave phenomena associated with novel acoustical properties.^1–4 Examples of novel acoustic metamaterial properties include zero or negative dynamic mass density⁵ and bulk modulus,⁶ chirality, and a more general material response that couples strain and momentum fields, as well as stress and velocity fields,^7–9 which is an acoustic analog to bianisotropy in electromagnetism.¹⁰ These novel properties open the door to extraordinary control of acoustic and elastic wave propagation including negative refraction,^11–13 superlensing,^14–16 cloaking,^17–19 extraordinary absorption,²⁰ and arbitrary control of the phase of waves reflected or transmitted from a sub-wavelength surface.^21–24 Early research on these topics often relied on sub-wavelength resonant metamaterial structures which limits their extraordinary behavior to narrow bands of frequency and whose performance can be degraded by the presence of absorption in their constituent materials.^1–4 Sub-wavelength structures that do not rely on resonances but lead to novel wave-bearing properties are therefore of significant interest for the future development of acoustic metamaterials and their use in practical applications. Some successful examples of non-resonant metamaterial structures to control acoustic and elastic waves include sub-wavelength elastic shells to produce spatially graded effective phase speed and impedance,^25,26 space-coiling structures to add an arbitrary phase to reflected or transmitted fields for field focusing and particle manipulation,^{21–24,27,28} and pentamode elastic materials for sound focusing and acoustic cloaking.^29–32 

One limitation of most metamaterial structures studied to date is that their dynamic material properties cannot be altered without changing the microstructure and thus the response of any given system cannot be changed after fabrication. However, there are many cases where the ability to change the dynamic effective material properties is highly desirable. For example, acoustic lenses cannot dynamically steer their focus and thus their use is more limited than electronically steered multi-element arrays. The ability to modulate the material properties of a medium in time and space has recently gained significant interest in the scientific community as it breaks parity-time symmetry and thus enables non-reciprocal acoustic and elastic wave phenomena.^33–36 One key challenge in useful material property modulation, or material tunability, is to design microstructures that enable sufficiently large changes in macroscopic response to be of interest in applications. Previous work on the topic of material property modulation or tuning has primarily focused on the use of piezoelectric elements controlled by electronic signals.^37,38 However, it is notable that the ability to change the material properties of a medium is inherently tied to its nonlinearity.⁴ Namely, if the material displays significant hardening or softening under finite deformation, then acoustic or elastic waves with infinitesimal amplitudes that propagate in the medium will do so at a speed that depends on the initial deformation. This behavior, often called small-on-large acoustic propagation, is of significant interest for detection and classification of damage using ultrasonic waves,^39,40 but has also been recently studied in the context of metamaterials.⁴¹ Methods developed by research in the latter category have also provided insights into static elastic non-reciprocity.⁴² Recent studies on material structures subjected to large deformation have focused on the absorption of large amplitude vibrations and impacts.^43–47 These types of structures capitalize on engineered structural instabilities to generate controlled and repeatable large deformation to collapse at specific force thresholds and capture the mechanical load, therefore providing impact isolation. When constrained, elastic instabilities can be represented as materials that demonstrate regimes of structural negative stiffness, where a reduction in stress levels leads to an increase in material deformation.^48–52 Of specific interest to the current work are efforts to design modified hexagonal honeycomb structures to absorb impacts via elastic buckling of curved beam-like structures in a repeating hexagonal lattice pattern known as a negative stiffness honeycomb (NSH).^43,44

Negative stiffness honeycombs are structured to generate a reversal of the usual force-displacement relationship when exposed to externally applied deformations. Specifically, these materials require a decreasing level of force when one imposes an increase in the displacement field.^43,44 This non-intuitive response is the result of structural instabilities at the unit cell level that result from geometric nonlinearity associated with large deflections of the constituent elastic elements in the structure. This behavior has been shown to be advantageous for the isolation from vibration and impacts.⁴⁵ In contrast to conventional honeycombs, negative stiffness honeycombs can recover to their initial shape after deformation, making the device reusable and therefore of interest for applications where repeated mechanical impacts may occur. While considerable research effort has been devoted to the study of these materials for impact isolation, there has been no investigation of these types of structures in regards to their elastic wave-bearing behavior. However, their strong geometric nonlinearity suggests that materials with similar microstructures are excellent candidates for tunable elastic metamaterials. This work therefore builds on these previous studies to investigate elastic wave motion in NSH and how this structure can be used as an acoustic metamaterial to control elastic wave propagation.

The finite element method is employed to solve the equations of motion in NSH due to the complex geometry and continuous nature of the NSH microstructure. The solution procedure involves two sequential steps. First, a macroscopic uni-axial deformation is applied to the NSH, resulting in increased stored potential energy and a new equilibrium configuration and geometry. The dynamic equations of motion are then linearized about the new equilibrium configuration and a Bloch wave analysis is used to determine the propagating modes along with their respective phase and group velocities as a function of direction. Bandgaps in materials that exhibit buckling can be tuned by pre-straining the material beyond the point of instability.^41,63 However, it is observed in the NSH that pre-strain below the point of instability is sufficient to provide significant tunability of bulk elastic wave propagation. The specific NSH geometry selected for the present analysis exhibits high levels of anisotropy, where the longitudinal mode is relatively insensitive to the pre-strain while other modes show very high sensitivity to the pre-strain. Due to this behavior, there exist partial bandgaps where only longitudinal wave motion exists in the elastic medium, therefore mimicking the behavior of pentamode elastic materials. In this case, the pentamode-like behavior can be tuned via a macroscopically imposed pre-strain. This work suggests that the pre-strained NSH may be one means of achieving tunable pentamode response with a level of anisotropy that can be controlled through the design of the microstructure.

The metamaterial unit cell under investigation is the modified version of the honeycomb design in Correa et al.⁴⁴ Like regular honeycombs, the NSH consists of an ordered configuration of prismatic cells. The NSH differs from conventional honeycomb structures in that NSH cells have been designed to permit a large, recoverable deformation and thus enable the absorption of energy impacts while returning to the original shape upon removal of the mechanical load. The recoverable deformation is the result of the unit cell geometry, which consists of curved beams that exhibit force-displacement behavior similar to that of bistable or snap-through structures.^43–45 A unit cell of the lattice under consideration is shown in Fig. 1, and the geometric parameters are shown in Fig. 2. The pre-curved beams are fabricated in the shape of a sinusoid with a spatial period L_x and an amplitude h_b, which mimics the shape of a beam buckled into its first Euler buckling mode.⁵³ The double pre-curved beam design as well as the vertical and horizontal elastic elements have been designed to mechanically constrain the system such that the system buckles into a pseudo-third buckling mode, which resembles a sinusoid with twice the spatial periodicity along the X₁-direction, when a specified external uni-axial pre-strain in the X₂ direction is attained.⁴⁴ In general, the elastic response of the NSH can be tuned through a few simple geometric parameters of the lattice, as shown in Fig. 2. The ratio $Q=hb/tb$ most directly affects the elastic stored energy in the unit cell; large ratios will yield a more pronounced snapping behavior.⁵³ In this work, a ratio of 4 is chosen to increase tunability, while it is noted that the unit cell will be bistable when Q is greater than 2.41.⁵³ Previous work has made use of the tunability of the lattice via these geometric features to produce a NSH that isolates vibration and mechanical impacts.⁴⁵ However, the strongly nonlinear elastic behavior of NSH materials suggests that linear wave propagation in these materials has a strong dependence on an externally imposed pre-strain. Further, because the large deformations are associated with buckling phenomena, the NSH can be fabricated from high-stiffness non-lossy materials without loss of the types of tunability described in this work. NSH are therefore of interest as candidate structures for elastic metamaterials that can control wave propagation by tailoring the dispersive nature of wave propagation through uni-axial pre-strain.

The NSH shown in Fig. 1 can be geometrically decomposed into a fundamental unit cell with a set of basis vectors in which the whole lattice can be reconstructed through translations of this fundamental unit cell along the basis vectors.⁵⁴ The unit cell chosen for this work is the shaded region in Figs. 1 and 2, which is spanned by the direct lattice vectors $a1,a2$⁠. Any spatially varying function in these materials must also obey the unit cell periodicity

where

and p₁, p₂ are integers that represent the unit cell location relative to the fundamental unit cell located at the position defined by $p1,p2=0$⁠. The direct lattice vectors are, in Cartesian coordinates

For periodic media, all unique propagating modes lie within the first Brillouin zone.⁵⁴ Therefore, the values of the Bloch wave vector $K$ are restricted to this region. The basis vectors in the reciprocal space are obtained through the following relationships:

where $z=a1×a2$⁠. From Eqs. (3) and (4), the reciprocal basis vectors are

The first Brillouin zone is then found graphically by creating a Wigner-Seitz cell with the above basis vectors.⁵⁴ Figure 3 shows the Brillouin zone overlaid on the NSH. Due to symmetry about the X₁ and X₂ axes, $K$ can be reduced to vectors that lie within the irreducible Brillouin zone, traced in black.

The total displacement field $ut$ in the NSH is assumed to be of the form

where $ud$ is the displacement due to the pre-strain, and $ua$ is the propagating acoustic wave, which is assumed to consist of small perturbations about the pre-strained displacement field (⁠$‖ua‖/‖ud‖≪1$⁠). This assumption allows for the computation of $ud$ and $ua$ in two sequential steps. First, $ud$ is found by solving the nonlinear static equation for a given pre-strain applied to the lattice. Then, $ua$ is computed by a linearized dynamic equation about the statically deformed state. A Lagrangian formulation is adopted for both modeling steps. Due to the periodicity of the lattice, only the displacements within a unit cell need to be computed for both problems, which greatly reduces the computational complexity.

The theory described in Secs. III A and III B is exact for pre-strain levels up to the point of instability.⁵⁵ Due to the geometric complexity, a finite element method is utilized for both the deformation and wave propagation studies. In this work, the open source code FEniCS is used to numerically perform the finite element algorithm.^56,57 The plane-strain approximation is used, reducing the computational domain to a two-dimensional surface.

Let $Ω⊂R2$ be the material domain of a representative unit cell in the undeformed configuration with coordinate X and boundary $Γ=Γp∪Γi$⁠, where Γ_p are the shared boundaries between Ω and the unit cell region in Fig. 1 and consists of the boundaries $ΓA+, ΓA−, ΓB+, ΓB−, ΓC+, ΓC−$ shown in Fig. 4, and Γ_i are the inner boundaries. Then, let $ϕ(X)$ be the deformation mapping, such that the displacement field is defined as

and the deformation gradient F defined as

An affine, uni-axial strain is applied to the NSH in the X₂ direction by application of a macroscopic deformation gradient $F¯$⁠, defined in matrix notation as⁶² 

where β is the macroscopic engineering pre-strain and the macroscopic strain $E¯=F¯TF¯−I=β2/2−β$ for $β<0$⁠. The microscopic deformation field is then given by the linear combination of the homogeneous deformation and a non-homogeneous deformation field $w(X)$⁵⁵ 

The microscopic deformation gradient is therefore given by

The micro- to macroscopic mapping of the deformation gradient leads to a constraint on $w$⁵⁵ 

where $n$ is the normal vector on the unit cell surface Γ_p. This constraint can be satisfied by requiring that $w$ satisfy periodic boundary conditions on the unit cell surfaces

where $r=A,B,C$ and the surfaces are specified in Fig. 4.

A variational formulation suited for a finite element algorithm is then needed to solve for $w$⁠, which is derived using the principle of stationary potential energy, written in the integral form as⁵⁸ 

where Π is the total potential energy, $W(F)$ is the strain energy density, and $t$ is the traction vector on the boundary. The interior boundaries Γ_i are assumed to be stress-free, and $t$ is assumed to be anti-periodic on Γ_p. Therefore, the boundary integral in Eq. (17) vanishes, leaving only the domain integral. The particular strain energy density function chosen in this work is the Saint Venant-Kirchhoff model, which takes into account the finite displacement but assumes small local strain values, and is written as⁵⁹ 

where λ and μ are the first and second Lamé parameters, respectively, E is the Green strain tensor, and I is the identity tensor. The use of this strain energy density function is acceptable for these lattices given the fact that the NSH have been designed to undergo elastic buckling and repeatable large deformation behavior of these systems has been experimentally demonstrated.^43,44

The equilibrium position is found by finding $w$ that minimizes Eq. (17). The minimum of Π can be found by defining a functional L that is equal to the directional derivative of Π with respect to an arbitrary vector $v$⁠, the so-called test vector in finite elements, and finding the root

where V is a suitable vector space with periodicity enforced. Equation (20) is discretized into a system of algebraic equations using FEniCS and Newton's method is used to find the root. This requires the Jacobian J, which is defined as

where $δw$ is a virtual displacement vector. Figure 5 shows the displacement in the lattice for $β=0.05$⁠, which was solved using the approach described above.

Once the displacement field $ud$ is determined from the variational principles discussed above, the acoustic displacement $ua$ is found by linearizing the nonlinear equations of motion about $ud$⁠. One way of accomplishing this is the use of Lyapunov's indirect method, which results in a linear wave equation with non-constant coefficients⁵⁸ 

where L is a spatially varying tensor known as the tangent modulus tensor, and ρ is the mass density. Due to the periodicity of the geometry and the coefficients L, the Floquet theorem is utilized to restrict the computational domain to the representative unit cell.⁴¹ The assumed form of the solution is thus given by

where $U(X)$ is a periodic function of the unit cell, and $eiK·X$ accounts for the phase change across each unit cell.⁴¹ The weak form of Eq. (22) is derived by the Hermitian inner product of the left- and right- hand sides with a test vector $v$⁠, written as

Application of Green's first identity then yields the following integral equation:

Modes propagating in the NSH are then determined by discretizing Eq. (26) into a system of equations, resulting in the matrix system

where the stiffness matrix A is a Hermitian matrix that is a function of the wavenumber $K$⁠, and the mass matrix B is positive-definite. Iso-frequency surfaces $ω(kx,ky)$ are generated by treating Eq. (27) as a generalized eigenvalue problem, where $ωn2$ and $Un$ are the eigenvalues and eigenvectors, respectively. For slowness and group velocity calculations, it is more useful to specify the frequency ω and calculate all wave vectors that exist as is shown in the following paragraph.

Let the wavenumber be expressed in polar form, $K=‖K‖[cos (θ) sin (θ)]T$⁠, where θ is the direction of propagation and $‖K‖$ is the wavenumber magnitude. Equation (27) can then be rewritten as

where $q=i‖K‖$⁠, and

The 0 and 1 subscripts in the above equations are to separate the k_x and k_y components of the sum in Eq. (26). Equation (28) is a quadratic eigenvalue problem for q, and can be solved by specifying ω and θ.

Once ω_n and $Un$ are obtained by solving Eq. (27), the group velocity vector $cg$⁠, with components

can be computed by taking the derivative of Eq. (22) with respect to k_x and k_y. The group velocity components are then given by⁶⁰ 

where $A,kx$ and $A,ky$ are the derivatives of the stiffness matrix with respect to the components of the wavenumber, which are given by

and

In this work, a NSH with geometric parameters given in Table I is investigated. The NSH is made of Nylon 11, with material properties shown in Table II.⁴⁴ Bloch wave solutions are found for four pre-strain values, $β=0.0, 0.0101,0.0207, 0.0252$⁠, which are illustrated in Fig. 6. These pre-strains were chosen such that the elastic wave propagation was stable for all wavenumbers. For higher pre-strain values, wave propagation may be unstable for one or more wavenumbers and thus is not considered. For each pre-strain, the tangent modulus tensor in Eq. (23) is computed, and the Bloch modes and natural frequencies are found by specifying the wave vector $K$ in the irreducible Brillouin zone and solving Eq. (27). To give insights into the physical vibration of the lattice, the modes are categorized into either a longitudinal, transverse, or higher order mode (H.O.M.). Once categorized, the modes can be tracked as a function of pre-strain, giving insights into the tunability of the phase speed of each mode, the allowable directions of wave propagation, and the direction of the propagation of energy carried by each mode.

The modes in the low wavenumber regime (near the Γ point) represent propagating plane waves whose wavelengths are very large when compared to the dimensions of the lattice. At the unit cell level, these modes correspond to the bulk motion of the unit cell moving either parallel to the wave front (longitudinal), perpendicular to the wavefront (transverse), or out-of-phase with neighboring unit cells (higher order modes). For a wave vector in this regime with a propagation angle of θ and a small but finite magnitude δk, these modes form a complete mass-orthogonal set $[UΓ,0(θ) UΓ,1(θ)⋯UΓ,N(θ)]$⁠, where the subscript Γ corresponds to modes near Γ, and are a function of the propagation direction, and N is the total number of degrees of freedom in the finite element model. In this work, $δk=0.01$ is chosen to ensure good separation of the longitudinal and transverse modes. Any mode $ϕ$ at a larger wavenumber in the first Brillouin zone in the direction θ can be decomposed into a sum of Γ-point modes

where c_i are the modal coefficients. These coefficients are determined using the mass-orthogonality of the Γ-point modes using the relation

For the results presented below, the classification of $ϕ$ is made using the Γ mode with the highest modal coefficient.

Figure 8(a) shows the dispersion branches on the outline of the irreducible Brillouin zone in Fig. 3 for the undeformed configuration. The mode on each branch is categorized using the modal decomposition technique discussed in Sec. IV A. From Fig. 8(a), it is evident that the mode categorization is not constant on a branch. This is due to mode coupling when two branches interact, which is common in an anisotropic media with high wavenumbers.⁶¹ Figures 8(b)–8(d) show the evolution of the dispersion curves as the pre-strain is increased. Overall, all branches decrease in frequency which is indicative of a strain-softening material. The result is that the effective phase speeds of each mode are reduced for increasing pre-strain. The most apparent changes are the decrease of the high-frequency higher order mode branches to a lower frequency range. The modes whose motion is parallel to the pre-strain loading direction are the most sensitive to the pre-strain, due to the decrease in the effective stiffness of the transversely loaded beam-like elements. This behavior is most noticeable with the transverse mode in the X₁ direction, whose displacement is parallel to the pre-strain loading. In addition, the longitudinal mode in the X₂ direction (whose displacement is also parallel to the pre-strain) is more sensitive to the pre-strain than the longitudinal mode in the X₁ direction. This can be seen by examining the ratio of the modal phase speed at a pre-strain level over the phase speed in the undeformed configuration at a set frequency. Figure 7(a) shows this in the X₁ direction at a frequency of 200 Hz, and Fig. 7(b) shows this with the same frequency in the X₂ direction. Indeed, the transverse wave phase speed in the X₁ direction changes by about 80% from the undeformed state to a pre-strain of 0.025, while the longitudinal wave only changes 4% over the same pre-strain range. The modal phase speeds in the X₂ direction exhibit a change of approximately 20% for the same range of pre-strains. Note that the transverse mode in the X₂ direction is also sensitive to the pre-strain, even though the motion is perpendicular to the loading direction. It is interesting to note that the longitudinal mode in the X₁ direction is the least sensitive of all modes. Also of significant interest is that for certain frequency ranges, only unimodal longitudinal wave motion is permitted in the lattice. For these frequency ranges, the lattice acts as a pentamode elastic material, or meta-fluid.³² The frequency range where this meta-fluid behavior occurs changes with pre-strain as discussed in Sec. IV C. The anisotropy of the NSH can be examined by investigating the slowness contours at a set frequency, which is calculated from the following relationships:

Figure 9(a) shows the slowness contours for the first three propagating modes in the first Brillouin zone for the undeformed configuration at a frequency of 1100 Hz. A horizontal dashed line at that frequency is provided for reference in Fig. 8. The undeformed configuration exhibits directional wave propagation at this frequency, namely the higher order branches are confined to propagate in directions nearly parallel to the X₂ direction while the transverse modes propagate at all angles except approximately $[32° 78°], [102° 148°], [212° 258°]$⁠, and $[282° 328°]$⁠. As the pre-strain is increased to $β=0.0101$⁠, shown in Fig. 9(b), the higher order modes can propagate in the X₁ direction, while the transverse mode can no longer propagate in the X₂ direction. Note that in Fig. 9, the longitudinal mode is very insensitive to the pre-strain. The normal vector to these slowness curves specifies the magnitude and direction of the group velocity, which is the direction of energy propagation, and is shown in Fig. 10. Figure 10(a) shows that pre-strain has little effect on the energy propagation of the longitudinal modes as the pre-strain is increased. However, the group velocity does become slower in the X₂ direction. Figure 10(b) is the transverse and higher order modes for the undeformed configuration. For this frequency, the shear waves carry energy in all directions, while higher order modes only carry energy near the X₂ direction. When a pre-strain of $β=0.0101$ is applied, the higher order modes can now carry energy in the X₁ direction, while the transverse mode can no longer carry energy in the X₂ direction.

The mode that is least sensitive to the pre-strain is the longitudinal mode. Because of this, the lattice can act as a meta-fluid for a specific range of frequencies and pre-strains in all directions. In the undeformed state, a meta-fluid region exists between approximately 1250 to 1750 Hz, shown in Fig. 8(a). As pre-strain is applied, the center frequency and bandwidth of this region decrease, as shown in Figs. 8(b)–8(d). The existence of meta-fluid behavior in all directions is confirmed by the slowness curve in Fig. 9(c). A pre-strained NSH can therefore mimic pentamode metamaterials, which are of interest for application of concepts from transformation acoustics such as acoustic cloaking¹⁷ and impedance-match gradient index lenses.³² In addition, the frequencies where this occurs as well as the degree of anisotropy can be tuned through the design of the unit cell and the pre-strain level, which will be reserved for future research.

In this work, a negative stiffness honeycomb is investigated as a potential elastic wave metamaterial that can be tuned through mechanical deformation. A small-on-large formulation is employed to model the total displacement in the NSH, which requires the solution of two sequential problems. First, the unit cell deformation due to an applied uni-axial pre-strain is computed. Once the new equilibrium configuration is found, a linearized wave equation about this equilibrium point is solved using a Bloch wave formulation, which yields the modes and frequencies of propagation. Due to the geometric complexity of the unit cell, a finite element method using the FEniCS code was utilized for both solution steps.

This lattice was shown to exhibit a large degree of tunability without requiring large pre-strains or strains that approach buckling instability configurations. The high-frequency higher order modes and the transverse mode in the X₁ direction were the most sensitive to the pre-strain. However, the longitudinal and transverse modes in the X₂ direction did exhibit a change of phase speeds of around 20%. In addition to the modal phase speeds, the direction of propagation for each mode was analyzed by investigating the slowness contours at a frequency of interest. The direction of energy propagation was inspected at the same frequency by calculating the group velocity. It was shown that the transverse and higher order modes were highly directional, and their direction of propagation can be tuned by changing the pre-strain amount. In contrast, the longitudinal modes propagated at all angles and were not sensitive to the applied pre-strain. One of the most interesting aspects of this NSH is its ability to act as a meta-fluid, supporting only longitudinal wave motion at all angles of propagation within a frequency band, which can be tuned with the pre-strain amount. The NSH can therefore mimic pentamode metamaterials in select frequency bands, which can be used in applications where an elastic material must be efficiently impedance-matched to an exterior fluid, such as an acoustic cloak or an elastic gradient index lens. The ability to tailor the geometric parameters and pre-strain amount, as well as its ability to recover to the initial configuration when unloaded, makes NSH a versatile metamaterial that can be used in numerous applications.

This material is based upon work supported by the National Science Foundation under Grant No. EFMA-1641078. The authors would also like to thank Sam Wallen for discussions on dispersion in periodic media.