One emerging research area within the fields of acoustic and elastic metamaterials involves designing subwavelength structures that display elastic instabilities in order to generate an effective medium response that is strongly nonlinear. To capture the overall frequency-dependent and dispersive macroscopic response of such heterogeneous media with subwavelength heterogeneities, a theoretical framework is developed that accounts for higher-order stiffnesses of a resonant, nonlinear inclusion that varies with a macroscopic pre-strain, and the inherent inertia associated with an inclusion embedded in a nearly incompressible elastic matrix material. Such a model can be used to study varying macroscopic material properties as a function of both frequency and pre-strain and the activation of such microscale instabilities due to an external, macroscopic loading, as demonstrated with a buckling metamaterial inclusion that is of interest due to its tunable and tailorable nature. The dynamic results obtained are consistent with similar static behavior reported in the literature for structures with elastic instabilities.

Acoustic and elastic metamaterials have attracted significant attention in recent years due to their ability to overcome limitations of conventional materials by exploiting engineered sub-wavelength structures.1 Traditionally, the fields of acoustic and elastic metamaterials focus on linear, periodic media with negative effective mass density or stiffness through design of subwavelength resonant structures, and the homogenization schemes developed to describe the effective properties. However, there is now an increasing interest in an effective material response that considers nonlinearity. Most studies of nonlinear acoustic metamaterials rely on spatial periodicity. Examples of periodic, one-dimensional nonlinear metamaterials include mass-spring chains where springs possess cubic nonlinearity,2 chains of bistable springs3 or beams,4 structured waveguides with periodic plates and side holes,5 and thin rods with finite internal mass-spring oscillators.6 Example nonlinear acoustic and elastic metamaterials in two dimensions include arrays of granular chains7 and deformable lattices with voided cavities,8 active components,9 or resonators containing coated masses, thin beams, and voids.10 Proposed applications include non-reciprocal devices,9,11,12 phononic switches and rectifiers,8,13–15 and tunable nonlinear lenses for filtering, focusing, or improved resolution of acoustic imaging devices.7,16,17

The acoustic metamaterials of interest generate the strong nonlinearity as a result of an engineered elastic compliance rather than by exploiting resonance phenomena. Some initial research on the nonlinear behavior of this type of metamaterial includes buckling beams in the form of negative stiffness honeycombs,18,19 woodpile structures,20 and compliant, three-dimensional (3D) elements known as buckliballs that are fabricated from soft polymers.21 Previous nonlinear acoustic metamaterial research has focused on large elastic deformation and dispersion analysis using Bloch–Floquet methods, which are limited to linear wave phenomena propagating through a structure subjected to large static deformation. However, heterogeneous media with randomly dispersed hyperelastic inclusions are potentially useful for exploiting traditional nonlinear wave phenomena such as phase conjugation, parametric amplification, Raman scattering, and enhanced energy dissipation.22 A model to study the homogenized frequency-dependent response of a metamaterial consisting of a dilute concentration of highly nonlinear subwavelength inclusions embedded in a viscoelastic material is developed here.

Example small-scale inclusions that gives rise to large microstructure-induced macroscopic nonlinearity are hollow encapsulated shells that exhibit buckling and are embedded in a host medium. These types of heterogeneous media have been studied for cases of fluid23,24 and soft elastic solid25–27 host media. Experiments have verified that large values of quadratic-order nonlinearity can be achieved as the microspheres buckle due to an external loading pressure in castor oil23 and nearly incompressible soft plastic.26 These results for buckling microspheres inspire work in buckling metamaterials. The advantage of buckling metamaterials is their potential to open up a larger design space for effective material properties, one that is more controllable and tailorable than the microspheres, and which may in turn push the limits of achievable nonlinearity.

The present work is influenced by the study of bubbly liquids and porous elastic media. Quasi-static homogenization methods for fluids and fluid-like media can be described using iso-pressure assumptions for both the effective medium and its constituents, such as the immiscible mixture law.28,29 Limits of the mixture law have been derived where the constituents possess small and approximately identical nonlinearity.30,31 For the case of a soft inclusion within a stiffer matrix, such as porous elastomers30 or bubbly liquids,28,32 it has been proven that only a small volume fraction is necessary to induce large nonlinearity of the effective medium. Therefore, it is anticipated that the strong nonlinearity of the metamaterial inclusions will exhibit strong macroscopic nonlinearity for dilute concentrations even when frequency dependence is considered.

Quasi-static assumptions are insufficient for describing an effective medium that contains viscoelastic constituents. Early work that took into account relaxation and nonlinearity associated with the viscoelastic elements led to wave equations that predict attenuation, dispersion, and harmonic generation in the resulting effective medium.33,34 In the present work, it is of interest to account for the effects of microscale inertia due to the dynamic oscillations of subwavelength inclusions, which can also influence the macroscopic linear and nonlinear properties. Similar dynamic homogenization methods have been applied to cases of air bubbles in water,35,36 as well as a soft elastomer matrix containing gas37 or voids.25 For air bubbles in water, viscous losses have also been taken into account.38 

These previously studied models, however, do not capture a more general nonlinear inclusion undergoing purely volumetric expansion, which is of interest to study small-scale metamaterial elements, nor do they account for the influence of microscale inertia of the heterogeneity, or a varying microscopic response due to a macroscopic loading condition. The present work develops a theoretical framework for a dynamic homogenization model that quantifies a dispersive effective medium comprised of compliant, nonlinear elements embedded in a nearly incompressible viscoelastic matrix material. The model is influenced by an effective medium model developed by Zabolotskaya and Soluyan for bubbly liquids.35,38 The present work not only extends the previously developed theory to account for an arbitrary spherical inclusion with strong material nonlinearity and non-negligible mass, but is also applied to an incremental response to account for finite deformations.

Expressions for incremental deformation theory are applied to the constitutive relationships of the inclusion and viscoelastic matrix in Sec. II. Two coupled equations are then derived in Sec. III to model the nonlinear wave propagation through a heterogeneous medium by following the analysis for gas bubbles in a fluid by Zabolotskaya and Soluyan.38 The frequency-dependent effective medium behavior is then obtained in Sec. IV by solving the coupled model equations at first- and second-order. Only second-harmonic generation is considered in this paper, but analyses for cubic nonlinearity are straightforward extensions of the present work.39 The low-frequency limit considered in Sec. V yields an evolution equation in the form of the Korteweg–deVries–Burgers equation. It is worth noting that the theory developed here is not only applicable to metamaterials. Instead, it is relevant for any microscale inclusions described by expanding the strain energy density to at most fourth-order in radial strain, which accounts for up to cubic nonlinearity in the dependence of stress on strain. For example, the present model could be relevant to buckling microspheres if an analytic constitutive relationship were given.

Also of interest in the present work is the frequency-dependent response and resulting harmonic generation, for example, metamaterial inclusion with designed elastic instabilities embedded in a nearly incompressible viscoelastic material proposed in Sec. VI. Such an element is not only tunable with an external activation, but tailorable with the geometric design and choice of material properties and may offer more controllable macroscopic material properties than its more conventional counterparts.

Consider a spherical inclusion that has a nonlinear constitutive relationship and is embedded in a viscoelastic matrix material. Of interest is the motion and deformation when accounting for both geometric and material nonlinearity of the inclusion and the matrix material. Finite deformations are therefore necessary to properly account for nonlinear effects. However, the constitutive behavior in nonlinear elasticity is more complicated than the linear elastic stress-strain relationship. To simplify the equations, incremental deformations from an applied pre-strain are introduced to describe both the inclusion and the matrix material. The small-on-large behavior at each increment represents a small perturbation about an imposed pre-strain.

The inclusion of interest is spherical, undergoing purely radial deformation due to uniform pressure imposed on the surface. Assume that the inclusion is at its global equilibrium configuration defined by radius R0 and pressure PI0. The inclusion is then deformed to a new reference state, defined as the pre-strain configuration through application of a static pressure PI1. The new configuration is defined by

(1)

which is a normalized change in radius from R0 to the pre-strained radius R1. In the present work, subscripts “1” will refer to the pre-strain states. Now, let the inclusion be further deformed a small amount to the total, instantaneous configuration, defined by pressure PI, radius R, and normalized change in radius

(2)

The difference between the total configuration and pre-strain is defined as the incremental deformation field. The incremental dimensionless radius is identical to the infinitesimal strain tensor for an object undergoing purely radial deformation. Therefore, let ε=ξξ1 be defined from Eqs. (1) and (2) and represent the strain increment from the reference state R1.

Since the deformation is radially symmetric, the stress tensor is only non-zero on the diagonal. For this particular case, the Cauchy stress tensor is related to the hydrostatic pressure with the relation σij=Pδij.40 The same definition holds when considering incremental deformations. Therefore, let the total pressure be expressed as a Taylor series expansion about the pre-strain. The incremental pressure pI may then be defined as

(3)

where the Taylor series coefficients, K1, K2, and K3, represent the linear, quadratically nonlinear, and cubically nonlinear stiffness moduli local to the pre-strain configuration and are defined as

(4)
(5)
(6)

The stiffness moduli K1, K2, and K3 all vary as a function of the applied pre-strain. Although Eq. (3) is obtained by assuming ε is a small quantity, the pressure-strain relationship has been expanded to cubic order because the proposed microscale inclusion possesses magnitudes of K2 and K3 that may be large enough to render the higher-order terms non-negligible.

It is assumed that the matrix is a nearly incompressible viscoelastic material. In this case, the shear stress in the matrix that interacts with the surface of the inclusion is accounted for through an effective pressure PM.41 The effective matrix pressure associated with the shear stress is formulated in terms of incremental deformations and evaluated at the surface of the inclusion as

(7)

where μM and AM refer to the local shear modulus and third-order elastic coefficient, which is a quadratic order term, respectively. These local moduli are defined explicitly based on Ref. 41 as

(8)
(9)

and PM1 denotes a pre-strain imposed on the matrix material. The following analysis is limited to μM0/KM01, where μM0 and KM0 are static values of the shear and bulk moduli, respectively, and AM0 is a static nonlinear elastic modulus at quadratic order.

A theoretical model is now derived to study nonlinear propagation of a pressure wave. It is assumed that the only source of acoustic nonlinearity is induced via the frequency-dependent oscillations of the inclusions, which accounts for material and geometric effects both of the inclusion and matrix material local to the inclusion surface. The change in volume of the inclusions due to an incident pressure wave represents a nonlinear source term within the otherwise linearized inhomogeneous wave equation. The derivation most closely follows that developed by Zabolotskaya and Soluyan as presented in Refs. 32 and 38 for bubbly liquids insofar as the main assumption is that the nonlinearity is dominated by the presence of the inclusions. As such, the inclusions appear as a distribution of volume (monopole) sources in an otherwise linear wave equation that is coupled to a nonlinear dynamical model of Rayleigh–Plesset form describing the response of the inclusions. While nonlinear propagation effects within the matrix itself may be worth including in subsequent work, they lie beyond the scope of the present work.

The volume of the entire effective medium V* is defined by the sum of the volumes of its constituents as V*=NVI+VM, where N is the number of inclusions, VI is the volume of the inclusion, and VM is the volume of the matrix material. Assume that the total volumes may be defined as the summation of a pre-strain, denoted with a subscript 1, and an increment. Divide the volume conservation equation by V*1 to obtain

(10)

The volume fraction at each pre-strain is defined as ϕ=NVI1/V*1, such that 1ϕ=VM1/V*1. The present analysis is limited to dilute volume fractions whereby ϕ1.

To obtain an approximation for the density of the effective medium, substitute the total and pre-strain densities for the effective medium volumes in Eq. (10). The total density may also be expressed as the summation of a pre-strain state and a small increment that represents an acoustic perturbation, where ρ=ρ1+ρ. The same substitution for density is made for volume of the matrix. Then, Eq. (10) may be expressed in terms of volume fraction, incremental effective densities of the effective medium and matrix, and incremental volume of the inclusion as

(11)

Only the linear terms in the binomial expansion of the density are needed to derive a linear wave equation with a source distribution characterized by the volume oscillations v of the inclusions. Binomially expand (1+ρM/ρM1)1 and (1+ρ*/ρ*1)1 up to linear order. The quadratic order terms in the expansions, (ρ*/ρ*1)2 and (ρM/ρM1)2 are neglected because their contributions are small relative to the strong material nonlinearity assumed for the inclusions. Such approximations have been reasonably made for bubbly liquids and should be valid for the current inclusions.32,38

Since the matrix material is nearly incompressible, ρM1ρM0ρM and the local stiffness moduli at each increment are identical to those at the global equilibrium. The density of the matrix is related to the pressure at linear order as ρM/ρM1p/ρMcMl2, where cMl2=[ KM0+(4/3)μM0 ]/ρM refers to the quasi-static, longitudinal sound speed of the surrounding matrix material. The nonlinearity in the equation of state of the matrix material is neglected because it is assumed to be much smaller than the nonlinearity of the inclusions. Common values of the parameter of nonlinearity at quadratic order, expressed as B/A, are on the order of 5–12 for fluids or fluid-like media.42 For example, B/A = 5 for water, and it is even smaller for air, for which B/A = 0.4. It is assumed that the strong nonlinearity of the inclusions of interest yields a value of B/A that is many orders of magnitude larger than that of the matrix material.39 In addition, it is also assumed that the matrix is much stiffer than the inclusions. For these assumptions, the contributions to the overall nonlinearity of the effective medium from the parameter of nonlinearity for the matrix material are negligible in comparison to that of the inclusion. Such an approximation is consistent with approximations made for bubbly liquids32,35 and is realized most simply for the immiscible mixture law for fluids in the quasi-static limit.29 It is therefore reasonable to assume a linearized equation of state for the matrix material.

Furthermore, it is assumed that the effective density is approximately that of the matrix material, such that ρ*1ρ*ρM. Although ρI is not necessarily small compared to ρM, reasonable values of ρI have little influence on the macroscopic density because ϕ1. With these simplifying approximations, Eq. (11) becomes

(12)

Since nonlinear convection terms are neglected in the derivation of Eq. (12), the linearized equations of mass and momentum conservation are combined to obtain a linear wave equation in the form of

(13)

Substitution of Eq. (12) into Eq. (13) to eliminate the density yields the following inhomogeneous wave equation:

(14)

Equation (14) is similar to that obtained in Refs. 32 and 38.

The change in volume of the inclusion is determined by an augmented form of the Rayleigh–Plesset equation including up to quadratic order nonlinearity. The current form describes the motion of an inclusion as defined in Sec. II and is derived following the formalism used in Ref. 41. This derivation appears in  Appendix A, where it is assumed that the surrounding matrix is a nearly incompressible viscoelastic material. The derived dynamics ignore interactions between inclusions, requiring the concentration to be dilute enough that each inclusion behaves independently. Given that ϕ1 is required to derive Eq. (14), it is reasonable that these interactions are ignored.

In terms of change in volume, the dynamic model equation up to quadratic order may be expressed as

(15)

Although Eq. (15) is of a similar form to the dynamical equation utilized for bubbly liquids, e.g., in Refs. 32 and 38, the coefficients defined in Eqs. (16)–(20) also account for the inertia and damping of the inclusion, contributions of a nearly incompressible elastic matrix, and the effects of the pre-strain.

The left-hand side of Eq. (15) describes a linear oscillator with undamped natural frequency (squared) given by

(16)

and the viscous loss factor

(17)

Both the undamped natural frequency ω1 and viscous loss factor δ are functions of pre-strain states and will vary as a function of the deformation through the variables ξ1, R1, ρI1, K1, and μM.

The right-hand side of Eq. (15) contains both the acoustic pressure acting on the inclusion and the quadratic nonlinearity associated with the oscillation of the inclusion. The coefficient of the acoustic pressure p in Eq. (15) is

(18)

The other two terms on the right-hand side of Eq. (15) account for quadratic nonlinearity. The coefficient b represents the nonlinearity in the constitutive relation acting on the surface of the inclusion, and includes both nonlinearity of the inclusion and the matrix material, where

(19)

In the linear, fluid limit for which μM = AM = 0, b is identical to the coefficient of nonlinearity β of the inclusion39,43 and it is similar to that obtained for bubbly liquids.38 Therefore, b can be considered to be an effective coefficient of nonlinearity that includes an increased stiffness and quadratic-order nonlinearity due to the effects of the shear stress of the matrix on the surface of the inclusion.

The remaining quadratic term is associated with inertia, the coefficient for which is defined as

(20)

The appearance of both ρI1 and ρM indicates that the inertial nonlinearity is associated with the moving mass of both the inclusion and the matrix material surrounding the inclusion, respectively. Although only quadratic nonlinearity is retained in the present analysis, accounting also for cubic order nonlinearity in Eq. (15) is relatively straightforward.39 

The coupled Eqs. (14) and (15) for p and v comprise the effective medium model. These equations can be solved together numerically or they can be solved iteratively to describe harmonic generation as in Sec. IV, where a specified form of the inclusion is considered. For acoustic waves with frequencies well below the natural frequency of the inclusion, a single nonlinear evolution equation in terms of the acoustic pressure can be obtained as in Sec. V.

Equations (14) and (15) can be solved iteratively to describe a wave at frequency ω that generates a secondary wave at frequency 2ω, following a standard expansion for the pressure and volume as done in Ref. 38 and subsequently reviewed in Ref. 32. The linear solution at frequency ω provides insight into the frequency-dependent attenuation and dispersion of the effective medium. The second-order solution at frequency 2ω describes the strength of the acoustic nonlinearity.

The pressure and volume are expressed in terms of the fundamental wave and its second harmonic as

(21)
(22)

It is assumed that the harmonic amplitudes pn and vn are both on the order of ϵn, where ϵ is typically defined as the acoustic Mach number and is a small parameter, such that ϵ1. Substitute Eqs. (21) and (22) into Eq. (14) and equate terms of order ϵ to obtain

(23)

Equation (15) then becomes

(24)

Now solve Eq. (24) for v1 and substitute that expression into Eq. (23) to obtain the following equation for the first-order pressure:

(25)

where

(26)

and

(27)

The quantity c̃1 in Eq. (25) is obtained by setting n = 1 in Eq. (26), but more generally c̃n is the complex phase speed at frequency . The real phase speed cn and attenuation coefficient αn at frequency are given by

(28)
(29)

Now consider second-harmonic generation. The wave equation at order ϵ2 is obtained by substituting Eqs. (21) and (22) into Eq. (14) to obtain

(30)

Likewise, Eq. (15) becomes

(31)

Combining Eqs. (24), (30), and (31) to eliminate v1 and v2 yields

(32)

where c̃2 is defined by Eq. (26) with n = 2 and

(33)

The magnitude of the complex function β̃(ω) provides a measure of the increase in acoustic nonlinearity introduced by the dynamic oscillations of the inclusions in connection with second-harmonic generation. There are two components of β̃, one involving the coefficient of constitutive nonlinearity b and the other involving the coefficient of inertial nonlinearity d multiplied by ω2. As noted by Zabolotskaya,36 these components may cancel if constitutive nonlinearity offsets inertial nonlinearity.

For propagation of a progressive plane wave in the +x direction with boundary conditions p1(x=0)=p0 and p2(x=0)=0 and source amplitude p0, the solution of Eq. (25) is

(34)

where k̃n=nω/c̃n is the complex wave number. Substitution of Eq. (34) into Eq. (32) yields a forced, inhomogeneous ordinary differential equation for p2, the solution of which is

(35)

The amplitude of the second harmonic varies with propagation distance from the source, frequency, and pre-strain. In the limit of weak dispersion, |k̃2|2|k̃1|, the amplitude of the second harmonic grows linearly with x. The regime of linear growth is the same as that identified for propagation through bubbly liquids.32,38

It may also be of interest to consider the low-frequency behavior that corresponds to the response of the effective medium at frequencies below the resonance frequency of the embedded inclusions. In the low-frequency limit, the higher-order inertial effects are much smaller than both the linear terms and the constitutive nonlinearity and may reasonably be discarded. Furthermore, the smallness of the quadratic terms and terms with time derivatives in Eq. (15) permit it to be rewritten recursively in the form of v=v(p) as follows:32 

(36)

Substitution of Eq. (36) into Eq. (14) then yields

(37)

Let the low-frequency sound speed be defined as

(38)

Recall that the parameters defining c* in Eq. (38) vary as a function of the pre-strain deformation and are therefore not constant. The term within the square brackets in Eq. (38) is the harmonic average of the compressibilities, defined as the inverse of stiffness, with ϕ1 and an effective stiffness of the inclusion (1+ξ1)[ K1+(4/3)μM ]. Therefore, c* represents the sound speed obtained by a simple mixture law. This approximation is applicable only to media subjected to iso-stress or iso-pressure, such as fluids with zero shear modulus or fluid-like media with shear a modulus much smaller than the bulk modulus.

A one-dimensional evolution equation for progressive plane waves is now derived. A slow scale with τ=tx/c* and x1=ϵx is introduced to reduce, at leading order, the wave equation in Eq. (37) to an evolution equation in the form of the Korteweg–deVries–Burgers equation:32 

(39)

where

(40)

denotes attenuation,

(41)

corresponds to dispersion, and

(42)

represents the quadratic acoustic nonlinearity in the low-frequency limit. As previously discussed, b represents an effective coefficient of nonlinearity that accounts for both the nonlinearity of the matrix associated with shear stress on the surface of the inclusion, and of the inclusion itself. This parameter is identical to βI in the fluid limit. Therefore, Eq. (42) is in the form that would be obtained for an immiscible mixture, with effective stiffness of the inclusion (1+ξ1)[ K1+(4/3)μM ], effective nonlinearity of the inclusion b, and negligible quasi-static nonlinearity of the matrix.29 

As an application of the theoretical framework developed above, the frequency-dependent effective medium behavior is obtained for hyperelastic inclusions embedded in a nearly incompressible viscoelastic medium that induces a nonlinear elastic response due to buckling of the inclusion structure. Sections VI A–VI C offer a qualitative understanding of the model capabilities and the behavior achievable for one example nonlinear inclusion with elastic instabilities. The inclusion geometry was chosen because the elastic deformation is known from previous studies to demonstrate strongly nonlinear behavior.44,45 However, the acoustic response due to propagating waves through a composite medium has not yet been explored, making the chosen design a prime candidate for the present work. The proposed metamaterial inclusion is appealing because the instabilities stem from the geometric design. Therefore, these instabilities are more tailorable and controllable than other buckling inclusions, such as buckling microspheres in a fluid23,24 or an elastic25–27 matrix material.

It should be emphasized that the theoretical model is valid not only for the chosen example inclusion, but for any class of inclusions, metamaterial or otherwise, that exhibit purely radial deformations and is reasonably described using Eqs. (1)–(6). It may also be used to solve the inverse problem to design, tailor, and optimize the inclusions for specific applications, which should be the focus of future work.

The chosen example metamaterial inclusion is shown in Fig. 1(a), where a vertical slice through the center of the sphere reveals a two-dimensional cut in Fig. 1(b). The inclusion is symmetric about both the cut plane and the highlighted symmetry lines. The inclusion features six sets of two orthogonal beams, where the thickness of each beam is a symmetric function of the arc length distance from its center. The beam is thinner in the center and on the ends that attach to the body of the inclusion in order to displace the beam more at those positions than at the thicker segments.44 The thickness modulated beam therefore behaves in the same manner as double beam elements that constrain the second buckling mode as described in Refs. 44 and 45. Elastic instabilities are induced by concentrating the imposed pressure at the inclusion-matrix interface at the center of the thickness modulated beam via the pressure transformer, as shown in Fig. 1. The inclusion is assumed to have an initial radius of R1 = 8.3 mm and is constructed out of nylon. The volume-averaged strain energy density is calculated using the commercial finite element software Comsol, where a displacement-controlled loading condition is imposed on the beam elements shown in Fig. 1 to induce the radial displacement of the inclusion.45 As is common practice in effective medium theory in general and in the field of metamaterials in particular, the complicated inclusion geometry depicted in Fig. 1 may be modeled as a homogeneous spherical inclusion with effective material properties subject to the restriction that the radius is much smaller than the wavelengths of interest. Therefore, the present work is also restricted to the long-wavelength limit. The strains, pressures, and incremental stiffness coefficients are then obtained numerically based on Eqs. (2)–(6) using techniques first described in Klatt and Haberman,45 which in turn describe the effective behavior of the homogenized spherical inclusion.

FIG. 1.

(Color online) An example metamaterial inclusion with a thickness modulated beam is displayed as a (a) full 3D sphere and (b) 3D sphere cut in half.

FIG. 1.

(Color online) An example metamaterial inclusion with a thickness modulated beam is displayed as a (a) full 3D sphere and (b) 3D sphere cut in half.

Close modal

The unstable nature of an element with buckled beams is briefly discussed here, but a more detailed description of the modeling of a similar geometry was provided by Klatt and Haberman.45 Although Klatt and Haberman considered double beam elements, the behavior is analogous to that of a thickness modulated beam.44 The pressure-versus-strain response is shown in Fig. 2(a), where the internal pressure of the inclusion PI appears in the solid line and the total pressure on the surface of the inclusion PI + PM is denoted with a dotted line. The stiffness versus strain appears in Fig. 2(b), where the local linear stiffness of the inclusion K1 is shown in the solid line and the effective stiffness when also accounting for contributions from the viscoelastic matrix K1+(4/3)μM is shown with a dotted line.

FIG. 2.

(a) Internal inclusion pressure PI in MPa (solid line) and total pressure on the surface of the inclusion PI+PM in MPa (dotted line) versus strain Ei and (b) local linear stiffness of the inclusion K1 (solid line) in MPa and total local linear stiffness of the inclusion and matrix K1+(4/3)μM in MPa (dotted line) versus strain EI. The three pre-strains of interest are denoted on the inclusion response with A (solid circle), B (solid triangle), and C (solid square) and on the inclusion plus matrix response as A (open circle), B (open triangle), and C (open square). The gray shaded region represents the microscale negative stiffness regime.

FIG. 2.

(a) Internal inclusion pressure PI in MPa (solid line) and total pressure on the surface of the inclusion PI+PM in MPa (dotted line) versus strain Ei and (b) local linear stiffness of the inclusion K1 (solid line) in MPa and total local linear stiffness of the inclusion and matrix K1+(4/3)μM in MPa (dotted line) versus strain EI. The three pre-strains of interest are denoted on the inclusion response with A (solid circle), B (solid triangle), and C (solid square) and on the inclusion plus matrix response as A (open circle), B (open triangle), and C (open square). The gray shaded region represents the microscale negative stiffness regime.

Close modal

The non-monotonic inclusion pressure versus strain is calculated numerically from the strain energy density and imposed displacement obtained from the finite element analysis via a displacement-controlled loading. However, a more realistic scenario for the inclusions is that of pressure loading. State A (solid circle) represents the global equilibrium state. The deformation of the inclusion itself, as opposed to one embedded in a matrix material, is determined by an external loading pressure that changes the internal pressure of the inclusion. As the externally imposed pressure increases, the inclusion compresses until it reaches the local maximum, denoted by the solid triangle at state B. This point denotes a strain state where the local linear stiffness is approximately zero, as evident by state B in Fig. 2(b).

The gray box in both Figs. 2(a) and 2(b) represents the microscale negative stiffness regime and corresponds to strain states where K10. When the inclusion is unconstrained, an additional increase in pressure would cause the inclusion to undergo a large and rapid reduction in radius, where it snaps to a higher pressure state rather than follow the solid black line. The behavior for unloading is similar. This so-called snap-through behavior is due to the elastic instabilities associated with the inclusion geometry. The inclusion configuration denoted by position C is unstable. Therefore, the inclusion cannot exist in that deformed state without a reaction force to counter the decreased pressure associated with an increased strain of the negative stiffness regime. For a fluid matrix where μM0 = 0 Pa, the total stiffness on the surface of the inclusion is only that of the inclusion.

In the context of wave propagation, these instabilities are often undesirable on the macroscale. Previous work has indicated that it is possible to obtain macroscopic stability for an unstable microscopic element if it is constrained in the microscale negative stiffness regime, which occurs when the shear modulus of the surrounding matrix material is sufficiently large relative to the stiffness of the unstable element.45–51 Macroscopic stability means that the curvature of the strain energy density function is positive for all macroscopic strains, which requires that the linear stiffness of the effective medium always be positive. This is achieved when K1+(4/3)μM>0. Therefore, the present example includes a matrix with a static shear modulus that is large enough in magnitude to obtain macroscopic stability of the heterogeneous medium at all deformation states, which requires that μM0>(3/4)K1/(1+ξ1) for all deformation states. If the shear modulus were not large enough, then the microscale negative stiffness would induce macroscale negative stiffness and the effective medium would exhibit instabilities manifested as snap-through phenomena similar to that of the inclusion. The present work not only confirms the static stability observed in other works45–51 but also offers a model to study the dynamic response at some stabilized static deformation. Macroscopic stability is therefore desired for the present example so that wave propagation can be considered for an inclusion constrained to any deformation state.

This model is also valid for a fluid matrix and will yield similar results for an inclusion with varying stiffness. However, the current case of constrained negative stiffness is considered more interesting because it offers the ability to tailor the macroscopic stiffness in a way that cannot be achieved with conventional, continuous materials and emphasizes the appeal of the proposed inclusion design.

Let the initial volume fraction of inclusion be only ϕ0 = 5 × 10−4. Although this volume fraction is very small, it is sufficient to induce large nonlinearity on the macroscale as discussed for soft inclusion within either a fluid matrix, such as for bubbly liquids,36 or within an elastic matrix with comparable quadratic nonlinearity to the inclusion.30,31 The matrix material is assumed to be a soft, rubber-like elastomer that is nearly incompressible,27,52,53 with the following material properties chosen to illustrate the behavior of interest: the bulk modulus is KM = 2.2 GPa, the density is ρM = 1000 kg/m3, and the shear modulus is μM0 = 1.5 MPa. To focus explicitly on the influence of the inclusion on the macroscopic behavior, AM0=0 Pa, such that the matrix has only linear material properties. Recall, however, that the effects of geometric nonlinearity have been included in the model and are still accounted for in the analysis presented here.

The effective pressure due to the shear stress of the matrix material is calculated numerically using Eqs. (7)–(9), where the strain increment ε in Eq. (7) is that of the inclusion. The total pressure acting on the surface of the inclusion PI + PM in Fig. 2(a) is now monotonic, which is indicative of macroscopic stability, and has an increased magnitude relative to PI. Thus, the value of μM0 is large enough to constrain the inclusions for all strain states, including those where K10. An imposed pressure results in loading and unloading that follows the dotted line, such that the inclusion is moving through the negative stiffness regime along the solid line. At pre-strain C, the shear stress of the matrix provides a sufficient reaction force to balance the decreased internal pressure of the inclusion for an increased external pressure, such that the beam elements of the inclusion remain in a deformed configuration that would be unstable in the absence of the restoring force provided by the deformed matrix. If the shear modulus were not sufficiently large, the reaction force would not constrain the inclusion, and the entire heterogeneous medium would exhibit snap-through phenomena associated with macroscopic instability. The effective local linear stiffness when accounting for the shear effects of the matrix becomes K1+(4/3)μM, as shown in the dotted line in Fig. 2(b), is always positive.

The frequency-dependent behavior is now calculated numerically based on the developed equations in Secs. III and IV for three pre-strain values depicted in Figs. 2(a) and 2(b), where an imposed pressure P0 constrains the inclusions to the states defined by A (circle), B (triangle), and C (square). Pre-strain A represents the global equilibrium. Although the pressure and strain at pre-strain A are identical, the stiffness differs in Fig. 2(b), where just the inclusion response is shown with the filled circle and the stiffness when accounting for the matrix appears with an open circle. Pre-strain B represents a point approximately at the zero linear stiffness state, but is a small, positive value. Pre-strain C denotes a state where the microscale inclusions are constrained to remain in a configuration that would be unstable without the presence of the matrix. This range is called the microscale negative stiffness regime.

In the present analysis, let the shear viscosity of the matrix be ηM=5Pas and the bulk viscosity of the inclusion be ζI=0Pas. Again, the matrix material property is within the range of values available for soft rubber-like materials in the literature for nearly incompressible media,52,53 and the viscosity of the inclusion is neglected to focus on other effects. From Eq. (17), it is apparent that the viscous damping due to the matrix material varies as a function of pre-strain even for a matrix displaying constant shear viscosity. Additionally, the undamped natural frequency from Eq. (16) also varies as a function of pre-strain. For ease of comparing the frequency-dependent properties of the effective medium at each pre-strain, it is advantageous to define a dimensionless frequency relative to one reference frequency. Let the undamped natural frequency at the zero pre-strain state, defined as ω0=ω1|ξ1=0, serve as the reference frequency for the discussion that follows. The frequency-dependent behavior for all pre-strain cases considered will be analyzed as a function of the dimensionless frequency ω/ω0.

In Fig. 3, the phase speed as defined in Eq. (26) is normalized by the quasi-static longitudinal sound speed of the matrix cMl and shown on a log scale. Well below the resonance frequency for all three cases, the phase speed c1 approaches the low-frequency sound speed c* defined by Eq. (38), which varies as a function of pre-strain. In the low-frequency limit, the effective medium stiffness is lower than that of the matrix material due to the presence of compliant inclusions oscillating in phase with the pressure wave. There is a decrease in stiffness at the pre-strain B (dotted-dashed line) relative to A (solid line), and thus the low-frequency phase speed for pre-strain B is less than A. Similarly, pre-strain C (dotted line) indicates a further decrease in overall stiffness and thus a lower phase speed.

FIG. 3.

(Color online) Normalized phase speed c1(ω)/cMl versus normalized frequency ω/ω0 for ϕ0=5×104 inclusions in an elastic matrix with μM0=1.5MPa at pre-strain A (solid line), pre-strain B (dotted-dashed line), and pre-strain C (dotted line).

FIG. 3.

(Color online) Normalized phase speed c1(ω)/cMl versus normalized frequency ω/ω0 for ϕ0=5×104 inclusions in an elastic matrix with μM0=1.5MPa at pre-strain A (solid line), pre-strain B (dotted-dashed line), and pre-strain C (dotted line).

Close modal

As the frequency increases, the phase speed decreases until the frequency equals the undamped natural frequency at each pre-strain, which corresponds to an antiresonance with minimum phase speed. As anticipated from Fig. 2(b), pre-strain C has the smallest natural frequency, then the natural frequency at pre-strain B, and pre-strain A case represents the largest natural frequency as exhibited in Fig. 3. For frequencies beyond the antiresonance, the phase speed variation becomes large, culminating in a maximum value at a frequency greater than ω/ω0 = 1. The large variation corresponds to frequencies in the range

(43)

The upper and lower bounds of the inequality in Eq. (43) are derived by considering the lossless case, where δ = 0. The lower frequency limit in Eq. (43) is obtained as c̃10 in Eq. (26), and the upper frequency bound corresponds to where c̃1 in Eq. (26). When damping is included, as in Fig. 3, the magnitude of c̃1 is finite, but may still become large within the bounds defined in Eq. (43). In the high-frequency limit beyond the upper bound of Eq. (43), c1(ω)cMl.

The frequency range defined in Eq. (43) represents a regime characterized by high attenuation. As the phase speed approaches infinity, which is the limit for an incompressible medium, the propagating waves attenuate more rapidly. Attenuation of the fundamental wave times the wavelength from Eq. (29) is shown in Fig. 4 for pre-strains A (solid line), B (dotted-dashed line), and C (dotted line). The wavelength λ1=c1/f is a function of the dispersive phase speed c1 from Eq. (28). In the frequency range where the phase speed variance is the largest, α1(ω)λ1 is non-zero and becomes large, but is very small outside of this frequency range. In Fig. 4, the peak magnitudes occur at the undamped natural frequency of each respective pre-strain. Pre-strain C possesses an attenuation coefficient with the largest width, indicating that the negative stiffness regime attenuates waves over the largest frequency range. Pre-strain B has the second largest bandwidth for high attenuation due to favorable damping properties of a system with zero linear stiffness state. The smallest bandwidth for high attenuation occurs for pre-strain A.

FIG. 4.

(Color online) Attenuation α1(ω)λ versus normalized frequency ω/ω0 for ϕ0=5×104 inclusions in an elastic matrix with μM0=1.5MPa at pre-strain A (solid line), pre-strain B (dotted-dashed line), and pre-strain C (dotted line).

FIG. 4.

(Color online) Attenuation α1(ω)λ versus normalized frequency ω/ω0 for ϕ0=5×104 inclusions in an elastic matrix with μM0=1.5MPa at pre-strain A (solid line), pre-strain B (dotted-dashed line), and pre-strain C (dotted line).

Close modal

The complex wavenumber k̃1 can be expressed as the summation of real and complex components, such that k̃1=ω/c1iα1. In the absence of loss in the matrix material or inclusion, the phase speed goes to zero at the natural frequency of each pre-strain, and the frequency range defined in Eq. (43) corresponds to no propagation. When ηM is small but non-zero, the phase speed is bounded but may still become large in magnitude, such that ω/c10 and k̃1iα1. In this case, waves are purely evanescent. However, if α1 is too large, the expansion for pressure and volume defined in Eqs. (21) and (22) are no longer valid since the assumptions |p2||p1| and |v2||v1| are no longer true. These regimes of high attenuation will be studied in future work.

Second-harmonic generation is used to illustrate the frequency and pre-strain dependence of the acoustic nonlinearity, which may become very large in magnitude. Let the following be defined as the normalized coefficient of quadratic nonlinearity:

(44)

where c̃2 and β̃ are the frequency and deformation dependent coefficients of nonlinearity and phase speed of the second harmonic defined in Eq. (26) (with n = 2) and Eq. (33), respectively. The low-frequency limits of c* and β*, as defined in Eqs. (38) and (42), respectively, are evaluated for the zero pre-strain case to provide a constant reference value when illustrating dependence on frequency and deformation.

The normalized coefficient of quadratic nonlinearity B(ω) appears as a function of frequency in Fig. 5. At all three pre-strains there exists two resonance peaks. The resonances are defined as the frequencies for which Δ1 = 0 and Δ2 = 0, where ωn appears in Eq. (27). For pre-strain A (solid line), one resonance occurs when the frequency equals the undamped natural frequency of the fundamental wave, such that ω/ω0 = 1. Another resonance peak occurs at a lower frequency, where ω/ω0 = 0.5, corresponding to the second harmonic frequency 2ω equaling the undamped natural frequency. The same qualitative behavior occurs for pre-strain B (dotted-dashed line) and pre-strain C (dotted line), but the frequencies of the resonances are shifted downward due to the decrease in undamped natural frequency at those pre-strains.

FIG. 5.

(Color online) Normalized coefficient of quadratic nonlinearity B(ω) versus normalized frequency ω/ω0 for ϕ0=5×104 inclusions in an elastic matrix with μM0=1.5MPa at pre-strain A (solid line), pre-strain B (dotted-dashed line), and pre-strain C (dotted line).

FIG. 5.

(Color online) Normalized coefficient of quadratic nonlinearity B(ω) versus normalized frequency ω/ω0 for ϕ0=5×104 inclusions in an elastic matrix with μM0=1.5MPa at pre-strain A (solid line), pre-strain B (dotted-dashed line), and pre-strain C (dotted line).

Close modal

For the case of gas bubbles in a liquid, analytic expressions for the second-harmonic amplitude were developed by Zabolotskaya and Soluyan38 for the cases of ω = ω0 and ω = 2ω0. Those expressions illustrate that there is growth of the second-harmonic amplitude with propagation distance up to a maximum value, and then a decrease in amplitude beyond that value as the gas bubbles absorb energy. Similar expressions are not derived in the present work, but the behavior for the present nonlinear inclusions will be qualitatively similar to that described for bubbly liquids and readers are referred to Ref. 38 for more information regarding such phenomena.

Although the same number of resonance peaks occurs at all three pre-strains, antiresonance exists only for pre-strain C. An antiresonance occurs at frequencies where β̃(ω)=0. Thus, the antiresonance occurs when the constitutive nonlinearity at quadratic order offsets the coefficient of inertial nonlinearity, i.e., for frequencies at which b=3(ω/ω1)2d in Eq. (33). When the coefficient of nonlinearity defined in Eq. (33) is negative, b is negative and only an imaginary frequency will result in b=3(ω/ω1)2d. Since an imaginary frequency is non-physical, no antiresonance exists when β̃(ω=0) is negative. The value of d is always positive, but no such restriction exists for b, which therefore determines the sign of β̃(0).

From Eqs. (4) and (5), the local quadratic stiffness moduli may be defined as a function of the local linear stiffness moduli through the following relation: K2K1/ξ. Thus, K2 can be approximated from the slope of the K1 versus EI curve shown in the solid line in Fig. 2(b). The slope of K1 is negative in the strain region between zero pre-stain and the local minimum of K1 and positive for negative strains greater in magnitude than the local minimum. Note that K2 is always negative for positive strains. For the specific example inclusion and matrix properties chosen, the sign of K2 determines the sign of both b from Eq. (19) and β* in Eq. (42). However, this will not always be the case; for example AM may offset K2 such that b is positive when K2 is negative. At pre-strain A and B, K2 is negative and thus β̃(0) is also negative. However, for the pre-strain C, both K2 and β̃(0) are positive and an antiresonance at ω=ω1b/3d is present.

In the low-frequency limit, the quadratic nonlinearity is dominated purely by the quasi-static nonlinearity of the inclusion and matrix defined by b in Eq. (19). The normalization in Eq. (44) is defined to ensure that the zero frequency limit of B(ω) at the global equilibrium is equal to 1, which is the case for pre-strain A. At pre-strain B, the normalized coefficient of quadratic nonlinearity at zero frequency is more than an order of magnitude larger than at pre-strain A, while pre-strain C is approximately 3 orders of magnitude larger than at pre-strain A in the low-frequency limit. Additionally, the peak amplitudes are largest for pre-strain C, second largest for pre-strain B, and smallest for pre-strain C. This illustrates that very large macroscopic acoustic nonlinearity can be induced due to the elastic instabilities on the microscale. Furthermore, the large range of variability possibly highlights the tunable nature of such metamaterial composites. In the high-frequency limit, where ω, the amplitude of the coefficient of nonlinearity decays as ω4 and β̃(ω)0. Therefore, beyond the undamped resonance at each pre-strain, the significance of the acoustic nonlinearity decreases as a function of increasing frequency.

Last, it is of interest to consider the amplitude of the fundamental and second harmonic as a function of distance from the source. The amplitude of the fundamental wave | p1 | (solid line), as calculated from Eq. (34), and second harmonic | p2 | (dotted line), obtained from Eq. (35), are shown in Pa over two propagation distances in Fig. 6. It is important that the behavior occurs over multiple wavelengths, and the lower x axis shows the dimensionless distance relative to the wavelength λ1. The upper x axis shows the distance from the source in meters. Although the dimensionless propagation distance x/λ1 is identical for Figs. 6(a)–6(c), the actual distance from the source differs for each case and the upper x axes are not identical for Figs. 6(a)–6(c). Propagation distance for which there is strong harmonic generation is relevant to the design and optimization of the inclusion.

FIG. 6.

(Color online) Propagation of | p1 | (solid line) and | p2 | (dotted line) over distance x/λ1 (lower x axis) and x in meters (upper x axis) from the source with p0=100 Pa at (a) pre-strain A with frequency ω/ω0=0.85 and nonlinearity |β̃|=1.7×102, (b) pre-strain B with frequency ω/ω0=0.34 and nonlinearity |β̃|=7.3×103, and pre-strain C with ω/ω0=0.23 and nonlinearity |β̃|=2.2×104. Note that the vertical range x in meters differs for each pre-strain.

FIG. 6.

(Color online) Propagation of | p1 | (solid line) and | p2 | (dotted line) over distance x/λ1 (lower x axis) and x in meters (upper x axis) from the source with p0=100 Pa at (a) pre-strain A with frequency ω/ω0=0.85 and nonlinearity |β̃|=1.7×102, (b) pre-strain B with frequency ω/ω0=0.34 and nonlinearity |β̃|=7.3×103, and pre-strain C with ω/ω0=0.23 and nonlinearity |β̃|=2.2×104. Note that the vertical range x in meters differs for each pre-strain.

Close modal

Since the phase speed, attenuation, and coefficient of nonlinearity vary with both frequency and strain, specific frequencies are chosen at each pre-strain of interest to explore the behavior of | p1 | and | p2 |. In Fig. 5, it was ascertained that the most prominent nonlinearity occurs at the undamped natural frequency of each pre-strain. However, high attenuation also occurs at the undamped natural frequency. Therefore, frequencies were chosen at ω/ω1 = 0.85 for each respective pre-strain to be far enough below resonance that attenuation does not completely dominate the response.

For pre-strain A, shown in Fig. 6(a), this frequency is ω/ω0  = 0.85 and the corresponding nonlinearity is |β̃|=1.7×102. The amplitude of the fundamental appears relatively constant as it propagates from the source. The magnitude of the second harmonic exhibits an oscillatory behavior with a peak amplitude more than 4 orders of magnitude smaller than | p1 |. The periodicity in | p2 | corresponds to the dispersion length 2π/| k22k1 |, where kn = ω/cn is the real part of the complex wave number. At distances less than one dispersion length, dispersion effects are negligible.

For pre-strain B, shown in Fig. 6(b), the chosen frequency of ω/ω1 = 0.85 now corresponds to ω/ω0 = 0.34 with |β̃|=7.3×103. The second harmonic still exhibits an oscillatory response over the period of one dispersion length, which is damped out due to the strong attenuation. The amplitude of the fundamental wave is also attenuated over the propagation distance shown. The peak amplitude of | p2 | is slightly more than 3 orders of magnitude smaller than | p1 |. The amplitude of the second harmonic at pre-strain B is larger than at pre-strain A, despite the higher attenuation, because |β̃| is an order of magnitude larger at pre-strain B than at A.

The lowest frequency considered is for pre-strain C, shown in Fig. 6(c), where ω/ω1 = 0.85 corresponds to ω/ω0 = 0.23 with nonlinearity |β̃|=2.2×104. This pre-strain also possesses the largest attenuation, and both the fundamental and second harmonic are attenuated as they propagate from the source more than at pre-strains A and B. The most attenuation is present at pre-strain C, as evident by the decaying amplitudes of both the fundamental and second harmonic. For | p2 |, the oscillatory response is damped out over the propagation distance shown more quickly than for pre-strain B. Pre-strain C also corresponds to the largest coefficient of nonlinearity, which is 1 order of magnitude larger than at pre-strain B and 2 orders of magnitude larger than at pre-strain A, and in turn, an increased amplitude of | p2 | in Fig. 6(c) relative to Figs. 6(a) and 6(b). At x/λ1=8,| p2 | is approximately 1 order of magnitude smaller than | p1 |. Thus, | p2 | from pre-strain C is closer in magnitude to | p1 | than for pre-strains A and B, indicating that pre-strain C generates the strongest second harmonic.

Even though |β̃| at pre-strain C is very large, | p2 || p1 | is still satisfied. For the case of strong dispersion, the phase speed of the second harmonic diverges from that of the fundamental, such that the quantity 1(c̃2/c̃1)2 in Eq. (35) may become large in magnitude. Thus, weak second harmonic generation relative to the fundamental wave may still be satisfied even for very large |β̃|. For weak dispersion, i.e., where |1(c̃2/c̃1)2|1,|β̃| is more restricted. With the current inclusion design, it may be possible to induce macroscopic nonlinearity that violates | p2 || p1 |, and future work should consider higher-order nonlinear effects.

The current work generalized the theoretical dispersion model developed by Zabolotskaya and Soluyan38 for bubbly liquids to account for spherical inclusions with a non-negligible moving mass within a viscoelastic matrix material subjected to macroscopic loading conditions. Both the micro- and macroscales are modeled assuming both material and geometric nonlinearity through incremental deformations. In addition to defining the frequency-dependent properties of an effective medium, including phase speed, attenuation, and coefficient of nonlinearity at quadratic order, the fundamental and second harmonic solutions for propagating waves are obtained. Furthermore, a low-frequency evolution equation in the form of the Korteweg–deVries–Burgers equation was derived.

The frequency-dependent wave propagation model also allows identification of three frequency regimes. At frequencies less than the undamped natural frequency of each pre-strain, the system is characterized by high nonlinearity. Harmonic generation is largest at these frequencies. The high-frequency range is where the inclusion motion is considered “frozen” and the effective properties are dominated by those of the matrix. An intermediate range between the low- and high-frequency regimes is where waves are strongly attenuated. The frequencies for each of these regimes depend on the pre-strain imposed on a low volume fraction of an example metamaterial inclusion. This result demonstrates the potential to vary the effective material parameters through mechanical loading of nonlinear acoustic metamaterials. Within the high-frequency limit, the metamaterial requirement that the microstructure be much smaller than the wavelengths of interest will be violated. Thus, the frequency-dependent wave propagation model also provides a means of bounding and verifying the long-wavelength limit, which is paramount to the study of metamaterials.

The present examples chosen reasonably neglected third-harmonic generation by requiring that | p2 || p1 |. By altering the inclusion design, matrix properties, or source amplitude, for example, cubic nonlinear effects and third-harmonic generation may become more prominent. Extensions to cubic order require accounting for cubic order constitutive and inertial nonlinearity coefficients in the volume formulation of a Rayleigh–Plesset type equation, which is a straightforward augmentation derived in the same way as the coefficient of nonlinearity at quadratic order.

The theoretical model presented in this paper is valid for any nonlinear inclusion that can be modeled following Sec. II and is not limited to the specific inclusion geometry with a thickness modulated beam used as the example here. However, the present model allows the frequency-dependent response for such a metamaterial inclusion to be studied. The frequency dependence of the phase speed, attenuation, and coefficient of nonlinearity all vary as a function of the applied pre-strain. The trends observed for the varying local linear stiffness of the specific inclusion of interest are consistent with behavior previously studied in the literature for a single structure or media that must be periodic. A key difference is that the dynamic response presented here is for a random arrangement of inclusions embedded in a matrix material.

As geometric and material parameters of the inclusion and matrix are altered, the qualitative behavior should remain the same, but the observed magnitudes will vary. For example, the magnitude of the attenuation may decrease, but the frequency range of large attenuation should still be broadest for an inclusion in the negative stiffness regime as compared to a zero linear stiffness state or a state where the constitutive curve is nearly linear. Thus, the theoretical model presented here offers a tool to study the macroscopic frequency- and loading-dependent properties induced by microscale nonlinear inclusions. The elements of interest here offer a means of creating a random effective medium that can yield enhanced attenuation and acoustic nonlinearity that can be tailored for specific applications.

This work was supported by the Applied Research Laboratories Chester M. McKinney Graduate Fellowship in Acoustics and a Multidisciplinary University Research Initiative from the Office of Naval Research (Grant No. N00014-13-1-0631).

Previous work has focused on the modeling gas bubbles in nearly incompressible elastic41 and viscoelastic54 media, as well as compliant objects in a nearly incompressible elastic material55 using Rayleigh–Plesset type equations. The present derivation combines the methodologies of these three papers, starting with Lagrange's equation,

(A1)

which is a function of kinetic energy T, potential energy U, dissipative function Ψ, and the instantaneous radius R. The overdot denotes derivatives with respect to time.

The matrix material is assumed to be nearly incompressible; therefore any change in the overall volume of the matrix is neglected. Thus, the kinetic energy is41 

(A2)

where ρM=ρM0 is the constant, equilibrium density and R denotes the instantaneous radius of the sphere. The kinetic energy for a compressible sphere undergoing purely radial deformation is

(A3)

There are three terms that contribute to the potential energy at the surface of a single inclusion: (i) the internal stored energy of the sphere, (ii) the effective potential energy due to the shear stress of the surrounding matrix, and (iii) the work done on the inclusion by an external force. It is assumed that potential energy is related to pressure such that41 

(A4)

The inclusion pressure is defined in Eq. (3). The effective shear pressure for incremental deformations is defined in Eq. (7) based on Ref. 41. Note that for a fluid, PM = 0. The third pressure is the work done by an external forcing function, and may contain time-harmonic and static components.

The dissipative function Ψ is derived from a dissipative function density ψ, related through the integral over volume. In the present work, losses are assumed to be linear. The dissipative function density for an isotropic solid with linear losses can be represented by56 

(A5)

where η and ζ are shear and bulk viscosity coefficients, respectfully, and ε̇ij is the infinitesimal strain rate, which neglects nonlinear dissipation terms. For the matrix, the only non-zero velocity is the radial component and only shear viscosity plays a role, where

(A6)

For the inclusion, only bulk viscosity is present due to the purely radial deformation, where

(A7)

It is assumed that both ηM and ζI are constants.

An extended Rayleigh–Plesset type equation is obtained by substituting Eqs. (A2)–(A4) and Eqs. (A6) and (A7) into Eq. (A1) and taking the necessary derivatives to obtain

(A8)

Now it is of interest to derive a volume formulation of Eq. (A8). First, let VI=VI1+v, where v is a small perturbation. Then, binomially expand the radius in powers of v up to quadratic order:

(A9)

Since the density ρI varies with deformation, it also must be expanded with respect to volume v, such that

(A10)

Note that the volume is expanded only up to linear order because it multiplies terms that are at least linear in order. Thus, for the product to be at most quadratic in order, higher-order terms for density are not necessary. Then, substitute Eqs. (A9) and (A10) into the left-hand side of Eq. (A8) and take the appropriate derivatives to obtain

(A11)

The loss terms are only retained to linear order because the expressions were derived assuming small strain rates. In terms of volume, the loss is quantified as

(A12)

For the pressure terms, ε must be expanded with respect to v, such that

(A13)

The combined internal pressure of the inclusion and effective shear pressure due to the matrix are defined in Eqs. (3) and (7) to obtain

(A14)

The example analysis only includes nonlinearity up to quadratic order, and thus K3 is not necessary for the present work. Extensions considering third-harmonic generation would require maintaining K3, as well as expanding all terms to cubic order in volume.39 Furthermore, assume that the external pressure can be expressed as the summation of an acoustic pressure wave and an external offset pressure, Pext=p+P0. To obtain the desired state P0=PI1+PM1.

From Eqs. (A11), (A12), and (A14), the dynamical equation for the volume up to quadratic order is

(A15)

where κ, ω1, δ, b, and d are defined in Eqs. (16)–(20).

1.
M. R.
Haberman
and
M. D.
Guild
, “
Acoustic metamaterials
,”
Phys. Today
69
(
6
),
42
48
(
2016
).
2.
K.
Manktelow
,
M. J.
Leamy
, and
M.
Ruzzene
, “
Multiple scales analysis of wave-wave interactions in a cubically nonlinear monoatomic chain
,”
Nonlinear Dynam.
63
,
193
203
(
2011
).
3.
N.
Nadkarni
,
C.
Darario
, and
D. M.
Kochmann
, “
Dynamics of periodic mechanical structures containing bistable elastic elements: From elastic to solitary wave propagation
,”
Phys. Rev. E
90
,
023204
(
2014
).
4.
J. R.
Raney
,
N.
Nadkarni
,
C.
Daraio
,
D. M.
Kochmann
,
J. A.
Lewis
, and
K.
Bertoldi
, “
Stable propagation of mechanical signal in soft media using stored elastic energy
,”
Proc. Sci. Acad. Sci. U.S.A.
113
(
35
),
9722
9727
(
2016
).
5.
L.
Quan
,
F.
Qian
,
X.
Liu
, and
X.
Gong
, “
A nonlinear acoustic metamaterial: Realization of a backwards-traveling second-harmonic sound wave
,”
J. Acoust. Soc. Am
139
(
6
),
3373
3385
(
2016
).
6.
R.
Khajehtourian
and
M. I.
Hussein
, “
Dispersion characteristics of a nonlinear elastic metamaterial
,”
AIP Adv.
4
,
124308
(
2014
).
7.
A.
Spadoni
and
C.
Daraio
, “
Generation and control of sound bullets with a nonlinear acoustic lens
,”
Proc. Natl. Acad. Sci. U.S.A.
107
(
16
),
7230
7234
(
2010
).
8.
K.
Bertoldi
and
M. C.
Boyce
, “
Wave propagation and instabilities in monolithic and periodically structured elastomeric materials undergoing large deformations
,”
Phys. Rev. B
78
,
184107
(
2008
).
9.
B.-I.
Popa
and
S. A.
Cummer
, “
Non-reciprocal and highly nonlinear active acoustic metamaterials
,”
Nat. Comm.
5
,
3398
(
2014
).
10.
P.
Wang
,
F.
Cadadei
,
S.
Shan
,
J. C.
Weaver
, and
K.
Bertoldi
, “
Harnessing buckling to design tunable locally resonant acoustic metamaterials
,”
Phys. Rev. Lett.
113
,
014301
(
2014
).
11.
B.
Liang
,
B.
Yuan
, and
J.-C.
Cheng
, “
Acoustic diode: Rectification of acoustic energy flux in one-dimensional systems
,”
Phys. Rev. Lett.
103
,
104301
(
2009
).
12.
X.
Guo
,
Z.
Lin
,
J.
Tu
,
B.
Liang
,
J.
Cheng
, and
D.
Zhang
, “
Modeling and optimization of an acoustic diode based on micro-bubble nonlinearity
,”
J. Acoust. Soc. Am.
133
(
2
),
1119
1125
(
2013
).
13.
N.
Boechler
,
G.
Theocharis
, and
C.
Daraio
, “
Bifurcation-based acoustic switching and rectification
,”
Nat. Mater.
10
,
665
668
(
2011
).
14.
F.
Li
,
P.
Anzel
,
J.
Yang
,
P. G.
Kevrekidis
, and
C.
Daraio
, “
Granular acoustic switches and logic elements
,”
Nat. Comm.
5
,
5311
(
2014
).
15.
J.-H.
Jang
,
C. H.
Koh
,
K.
Bertoldi
,
M. C.
Boyce
, and
E. L.
Thomas
, “
Combining pattern instability and shape-memory hysteresis for phonic switching
,”
Nano Lett.
9
(
5
),
2113
2119
(
2009
).
16.
C. M.
Donahue
,
P. W. J.
Anzel
,
L.
Bonanomi
,
T. A.
Keller
, and
C.
Daraio
, “
Experimental realization of a nonlinear acoustic lens with a tunable focus
,”
App. Phys. Lett.
104
,
014103
(
2014
).
17.
B. M.
Goldsberry
and
M. R.
Haberman
, “
Negative stiffness honeycombs as tunable elastic metamaterials
,”
J. Appl. Phys.
123
,
091711
(
2018
).
18.
D. M.
Correa
,
T.
Klatt
,
S.
Cortes
,
M.
Haberman
,
D.
Kovar
, and
C.
Seepersad
, “
Negative stiffness honeycombs for recoverable shock isolation
,”
Rapid Prototyping J.
21
(
2
),
193
200
(
2015
).
19.
S.
Shan
,
S. H.
Kang
,
J. R.
Raney
,
P.
Wang
,
L.
Fang
,
F.
Candido
,
J. A.
Lewis
, and
K.
Bertoldi
, “
Multistable architected materials for trapping elastic strain energy
,”
Adv. Mater.
27
,
4296
4301
(
2015
).
20.
E. B.
Duoss
,
T. H.
Weisgraber
,
K.
Hearon
,
C.
Zhu
,
W.
Small
, IV
,
T. R.
Metz
,
J. J.
Vericella
,
H. D.
Barth
,
J. D.
Kuntz
,
R. S.
Maxwell
,
C. M.
Spadaccini
, and
T. S.
Wilson
, “
Three-dimensional printing of elastomeric, cellular architectures with negative stiffness
,”
Adv. Fund. Mater.
24
,
4905
4913
(
2014
).
21.
J.
Shim
,
C.
Perdigou
,
K. R.
Chen
,
K.
Bertoldi
, and
P. M.
Reis
, “
Buckling-induced encapsulation of structured elastic shells under pressure
,”
Proc. Natl. Acad. Sci. U.S.A.
109
(
16
),
5978
5983
(
2012
).
22.
F. V.
Bunkin
,
Y. A.
Kravtsov
, and
G. A.
Lyakhov
, “
Acoustic analogues of nonlinear-optics phenomena
,”
Sov. Phys. Usp.
29
,
607
619
(
1986
).
23.
D. H.
Trivett
,
H.
Pincon
, and
P. H.
Rogers
, “
Investigation of a three-phase medium with a negative parameter of nonlinearity
,”
J. Acoust. Soc. Am.
119
(
6
),
3610
3617
(
2006
).
24.
P.
Marmottant
,
S.
van der Meer
,
M.
Emmer
,
M.
Versluis
,
N.
de Jong
,
S.
Hilgenfeldt
, and
D.
Lohse
, “
A model for large amplitude oscillations of coated bubbles accounting for buckling and rupture
,”
J. Acoust. Soc. Am.
118
(
6
),
3499
3505
(
2005
).
25.
L. A.
Ostrovsky
, “
Wave processes in media with strong acoustic nonlinearity
,”
J. Acoust. Soc. Am.
90
(
6
),
3332
3337
(
1991
).
26.
V. Y.
Zaitsev
,
A.
Dyshkin
,
E.
Pasternak
, and
L.
Matveev
, “
Microstructure-induced giant elastic nonlinearity of threshold origin: Mechanism and experimental demonstration
,”
Europhys. Lett.
86
,
44005
(
2009
).
27.
R.
de Pascalis
,
I. D.
Abrahams
, and
W. J.
Parnell
, “
Predicting the pressure-volume curve of an elastic microsphere composite
,”
J. Mech. Phys. Solids
61
,
1106
1123
(
2013
).
28.
E. C.
Everbach
, “
Tissue composition determination via measurement of the acoustic nonlinearity parameter
,” Ph.D. thesis,
Yale University
,
1989
.
29.
E. C.
Everbach
,
Z.
Zhu
,
P.
Jiang
,
B. T.
Chu
, and
R. E.
Apfel
, “
A corrected mixture law for B/A
,”
J. Acoust. Soc. Am.
89
(
1
),
446
447
(
1991
).
30.
I. Y.
Belyaeva
and
V. Y.
Zaitsev
, “
Nonlinear elastic properties of microinhomogeneous hierarchically structured media
,”
Acoust. Phys.
43
(
5
),
510
515
(
1997
).
31.
I. Y.
Belyaeva
and
V. Y.
Zaitsev
, “
The limiting value of the parameter of elastic nonlinearity in structurally inhomogeneous media
,”
Acoust. Phys.
44
(
6
),
635
640
(
1998
).
32.
M. F.
Hamilton
,
Yu. A.
Il'inskii
, and
E. A.
Zabolotskaya
, “
Dispersion
,” in
Nonlinear Acoustics
, edited by
M. F.
Hamilton
and
D. T.
Blackstock
(
Acoustical Society of America
,
Melville, NY
,
2008
), Chap. 5, pp.
167
175
.
33.
V. Y.
Zaitsev
,
V. E.
Nazarov
, and
I. Y.
Belyaeva
, “
The equation of state of a microinhomogeneous medium and the frequency dependence of its elastic nonlinearity
,”
Acoust. Phys.
47
(
2
),
178
183
(
2001
).
34.
V. E.
Nazarov
,
V. Y.
Zaitsev
, and
I. Y.
Belyaeva
, “
Nonlinear transformation of acoustic waves in microinhomogeneous media with relaxation
,”
Acta Acust. united Acust.
88
(
1
),
40
49
(
2002
).
35.
E. A.
Zabolotskaya
and
S. I.
Soluyan
, “
Emission of harmonic and combination frequency waves by air bubbles
,”
Sov. Phys. Acoust.
18
(
3
),
396
398
(
1973
).
36.
E. A.
Zabolotskaya
, “
Acoustic second-harmonic generation in a liquid containing uniformly distributed air bubbles
,”
Sov. Phys. Acoust.
21
(
6
),
569
571
(
1976
).
37.
L. A.
Ostrovskii
, “
Nonlinear acoustics of slightly compressible porous media
,”
Sov. Phys Acoust.
34
(
5
),
523
526
(
1988
).
38.
E. A.
Zabolotskaya
and
S. I.
Soluyan
, “
Nonlinear wave propagation in a liquid containing uniformly distributed air bubbles
,”
Sov. Phys. Acoust.
19
(
5
),
442
444
(
1974
).
39.
S. G.
Konarski
, “
Tunable, nonlinear acoustic metamaterials due to subwavelength structural instabilities
,” Ph.D. thesis,
The University of Texas at Austin
, Texas,
2017
.
40.
Y. C.
Fung
and
P.
Tong
,
Classical and Computational Solid Mechanics
(
World Scientific Publishing Co. Pte. Ltd.
,
Singapore
,
2001
), pp.
118
119
, 139–140, 421.
41.
S. Y.
Emelianov
,
M. F.
Hamilton
,
Yu. A.
Ilinskii
, and
E. A.
Zabolotskaya
, “
Nonlinear dynamics of a gas bubble in an incompressible elastic medium
,”
J. Acoust. Soc. Am.
115
(
2
),
581
588
(
2004
).
42.
R. T.
Beyer
, “
The parameter B/A
,” in
Nonlinear Acoustics
, edited by
M. F.
Hamilton
and
D. T.
Blackstock
(
Acoustical Society of America
,
Melville, NY
,
2008
), Chap. 2, pp.
25
37
.
43.
R. N.
Thurston
and
M. J.
Shapiro
, “
Interpretation of ultrasonic experiments on finite-amplitude waves
,”
J. Acoust. Soc. Am.
41
(
4
),
1112
1125
(
1967
).
44.
J.
Qiu
,
J. H.
Lang
, and
A. H.
Slocum
, “
A curved-beam bistable mechanism
,”
J. Microelectromech. Sys.
13
(
2
),
137
146
(
2004
).
45.
T.
Klatt
and
M. R.
Haberman
, “
A nonlinear negative stiffness metamaterial unit cell and small-on-large multiscale material model
,”
J. Appl. Phys.
114
,
033503
(
2013
).
46.
R. S.
Lakes
and
W. J.
Drugan
, “
Dramatically stiffer elastic composite material due to a negative stiffness phase?
,”
J. Mech. Phys. Solids
50
,
979
1009
(
2002
).
47.
W. J.
Drugan
, “
Elastic composite materials having a negative stiffness phase can be stable
,”
Phys. Rev. Lett.
98
,
055502
(
2007
).
48.
D. M.
Kochmann
and
W. J.
Drugan
, “
Analytical stability conditions for elastic composite materials with a non-positive-definite phase
,”
Proc. R. Soc. A
468
(
2144
),
2230
2254
(
2012
).
49.
C. S.
Wojnar
and
D. M.
Kochmann
, “
A negative-stiffness phase in elastic composites can produce stable extreme effective dynamics but not static stiffness
,”
Philos. Mag.
94
(
6
),
532
555
(
2013
).
50.
D. L.
Platus
, “
Negative-stiffness-mechanism vibration isolation system
,”
Proc. SPIE
1619
,
44
54
(
1991
).
51.
L. B.
Kashdan
,
C. C.
Seepersad
,
M. R.
Haberman
, and
P. S.
Wilson
, “
Design, fabrication, and evaluation of negative stiffness elements using sls
,”
Rapid Prototyping J.
18
(
3
),
194
200
(
2012
).
52.
A. M.
Baird
,
F. H.
Kerr
, and
D. J.
Townsend
, “
Wave propagation in a viscoelastic medium containing fluid-filled microspheres
,”
J. Acoust. Soc. Am.
105
(
3
),
1527
1538
(
1999
).
53.
C. C.
Church
, “
The effects of an elastic solid surface layer on the radial pulsations of gas bubbles
,”
J. Acoust. Soc. Am
97
(
3
),
1510
1521
(
1994
).
54.
E. A.
Zabolotskaya
,
Yu. A.
Ilinskii
,
G. D.
Meegan
, and
M. F.
Hamilton
, “
Modification of the equation for gas bubble dynamics in a soft elastic medium
,”
J. Acoust. Soc. Am.
118
(
4
),
2173
2181
(
2005
).
55.
E. A.
Zabolotskaya
,
Yu. A.
Ilinskii
, and
M. F.
Hamilton
, “
Weakly nonlinear oscillations of a compliant object buried in soil
,”
J. Acoust. Soc. Am.
125
(
4
),
2035
2040
(
2008
).
56.
E. A.
Zabolotskaya
,
M. F.
Hamilton
,
Yu. A.
Ilinskii
, and
G. D.
Meegan
, “
Modeling of nonlinear shear waves in soft solids
,”
J. Acoust. Soc. Am.
116
(
5
),
2807
2813
(
2004
).