This paper derives and demonstrates a one-dimensional acoustic metamaterial homogenization method. The homogenization method uses a multiple-scales approximation with Hamilton's principle, a weak-form representation of the dynamic equation. While the multiple-scales approximation makes the predicted effective material properties of this method inexact, the method is shown to be highly versatile. Analytical and numerical examples are given showing the ability of the homogenization method to account for viscosity and finite-amplitude effects.

Metamaterial research, or the effort to create synthetic composite materials with specific, often exotic, effective properties, relies upon homogenization methods. Homogenization methods are ways to replace inhomogeneous systems with an effective homogeneous system for analysis purposes. Many homogenization methods for acoustic and elastic metamaterials have been proposed. These methods range from analytical schemes such as scattering-based methods1,2 and statistical mechanics (see, for example, Ref. 3) and numerical methods such as self-consistent models,4 to experimental methods.5,6 The majority of the theoretical methods describe micromechanical behavior using a strong-form, or differential, representation rather than a weak-form, or integral, representation. There are many reasons for using a strong-form representation. The physical interpretation of strong forms is relatively easy to understand, and most existing analytical descriptions of fluid and solid mechanics are presented in differential form. On the other hand, weak-form representations, such as Hamilton's principle, have traditionally been used in high level derivations such as in the Lagrangian and Hamiltonian formulation of mechanics and, more recently, in finite element-based numerical methods.

The fact that integration inherently considers finite-sized domains makes using a weak-form representation of the microstructure attractive from a homogenization standpoint. In fact, a statistical mechanics analysis of a perfect gas proceeds from Hamiltonian mechanics. A more recent example is the work by Andia et al.,7 where a Lagrangian-mechanical approach was developed with the intention to incorporate molecular dynamics into their homogenization scheme. Their work focused primarily on periodic systems, and obtains effective properties by taking the variation of the microscopic Lagrangian with respect to the averaged (macroscopic) displacement.

The purpose of this paper is to present a new, multi-scale homogenization method based on Hamilton's principle for mechanical media. Hamilton's principle is, as mentioned above, a weak-form representation and is equivalent to Newtonian mechanics for continuous media. An important feature of Hamilton's principle is the ease with which constraints, such as constitutive equations and complicated geometries, may be imposed.

The outline of this paper proceeds as follows. Section II presents a basic introduction to Hamilton's principle and then derives the new homogenization method. Examples that illustrate aspects of the new homogenization method are then given in Sec. III. The examples include accounting for viscous damping and considering homogenization of material properties associated with finite-amplitude propagation. Finally, Sec. IV summarizes the paper.

Newton's second law for one-dimensional acoustic systems may be cast in weak form by Hamilton's principle, which takes the form8,9

(1)

where δ denotes the variation, S is the cross-sectional area (assumed to be constant throughout the analyses below), Ekin is the kinetic energy density of the system, Epot is the potential energy density of the system, F is the sum of all body force densities (which may be a function of space, time, and acoustic fields), and u is the particle displacement. The difference of the kinetic and potential energy densities is also called the Lagrangian density. Hamilton's principle is an optimization problem, and so additional constraints may be added via a Lagrange multiplier. For example, incorporating the constraint C = 0, where C is some function of the acoustic variables, may be done by adding the term SμC to the Lagrangian density, yielding a modified equation,

(2)

Here, the yet-to-be determined function μ acts as a Lagrange multiplier.

For homogeneous and linear acoustic systems, the kinetic and potential energy densities may take the form, respectively, of

(3)
(4)

where ρ is the mass density, β is the linear compressibility, v=u/t is the particle velocity, and p is the acoustic pressure. A relevant constraint is the constitutive equation

(5)

Using these definitions and assuming no external body forces (i.e., F = 0; this restraint will be relaxed below), Eq. (2) becomes

(6)

Carrying out the variation then yields

(7)

The variations of the particle velocity and pressure are assumed to be smooth functions of space. Then integrating by parts to isolate the variations of the acoustic pressure and particle velocity yields

(8)

Since variations are arbitrary the integrand of the left-hand side of the above equation yields three independent equations that must be simultaneously satisfied in the bulk medium:

(9)

The last of these equations is just the constraint C = 0, and the first two may be used to eliminate the Lagrange multiplier μ, yielding

(10)

which is just Newton's second law recovered.

The initial and final terms on the right-hand side of Eq. (8) may be neglected by assuming the pressure is zero at the possibly remote times t0 and t1. The boundary term may be re-written as

(11)

As is shown below, this form of the boundary term shows that u and p must be independently continuous at interfaces. Consider an integral over x(x0,x1). This integral may be arbitrarily broken into two integrals, one over x(x0,x) and one over x(x,x1), where x0<x<x1. Each of these integrals lead to boundary conditions of the same form as Eq. (11), leading to the corresponding equation

(12)

where the subscript < denotes the quantity approaching x with x<x and the subscript > denotes the quantity approaching x with x>x. Since the final term on the right is identical to the boundary term without splitting the integral, we find that the remaining terms on the right must be equal to zero. This can only be the case if p<δu<=p>δu>. Since the variation is arbitrary, we may therefore conclude that p<=p> and δu<=δu>, or, since the variation is also assumed to be smooth, u<=u>. The two integrals may be treated as separated domains that might have had independent material properties. Thus, the above result implies that the acoustic pressure and the particle displacement are continuous at arbitrary interfaces.

Hamilton's principle as described above is a global description of the system, or every part of the system must be accounted for simultaneously. However, the Lagrangian density and the constraint equation are written in terms of local variables. Thus, Hamilton's principle provides a way to transition from locally valid approximations to a globally valid approximation. Section II B presents a multiple scales-based analysis providing locally valid approximations of the acoustic variables, and Sec. II C explains the globally valid homogenization.

Consider a static, layered, inhomogeneous, one-dimensional, lossless acoustic system with a representative volume element length of L. Each layer is assumed to be a homogeneous acoustic/elastic material. By defining L we implicitly assume that the statistics of the inhomogeneities are stationary and ergodic, or that in any sub-domain of length L of the system the relevant statistics are, approximately, the same. (Note that choosing L for a given system is a non-trivial problem; see Ref. 10 for a discussion.) Locally, the behavior of the system may be represented by the constitutive equation and the equation of conservation of momentum, which may be written (in Lagrangian coordinates), respectively, as

(13)
(14)

where

(15)

is the strain and B is a coefficient of quadratic nonlinearity in the constitutive relation. The ellipsis at the end of the constitutive equation denotes that higher-order terms in the strain may also appear. To determine the relative importance of each of these terms, define the characteristic amplitudes of the pressure and velocity p0 and v0, respectively, and the smallest characteristic time scale of the acoustic signal T. These characteristic quantities are defined such that Πp/p0 and Λv/v0 are O(1) and derivatives of acoustic variables with respect to time are O(T). Defining σx/L and τt/T, Eqs. (13)–(15) may then be written as

(16)
(17)
(18)

where ν/τΛ. Let the characteristic acoustic Mach number be ϵv0/cmin, where cmin is the smallest small-signal sound speed present in the system. Assuming weak nonlinearity such that terms O(ϵ2) relative to the largest order term may be neglected, we may write

(19)
(20)

where ηL/cminT,Mβp0cmin/v0,Rρcminv0/p0,N1+B.

For η1 and η2ϵ=O(1) (long-wavelength limit with relatively weak nonlinearity) a multiple-scales analysis may be used to perform an asymptotic homogenization. Define the inhomogeneity-scale position σ0σ and the nth macroscopic position σnσ/ηn such that

(21)

Assume that any acoustic field X (e.g., Λ or Π) may be written as a convergent power series in η as X=X0+ηX1+η2X2+. Substituting these series forms into Eqs. (19) and (20) and collecting powers of η yields

(22)
(23)

The coefficients of each order of η may be treated as independent terms leading to separate equations. The O(η0) equations are

(24)

and the O(η1) equations are

(25)

Note that the nonlinear terms do not appear in the O(η1) equations due to the O(η0) equations.

Since Π0 and Λ0 are independent of σ0 the O(η1) acoustic variables will contain secular terms (i.e., unbounded in σ0) if the right-hand sides of Eqs. (25) are not identically equal to zero, yielding

(26)

These equations may be combined yielding the linear wave equations accurate to O(η2) (see Ref. 11):

(27)

Note that while the effective wave speed is a function of σ0 and is inhomogeneous, the wave equation is homogeneous in σ1. Any solution which is smooth in τ (assumed to be true since the smallest characteristic time scale is T > 0) is also smooth in σ1, and so a convergent Taylor series may be written as

(28)

Recalling that σ1=σ/η then shows that the first term is O(η0), the second term is O(η1), the third term is O(η2), and so on. However, the solution X=X0 as derived above is only valid to O(η2), and so we may write

(29)

Thus the acoustic variables may be considered completely independent of the microstructure to O(η1), and vary only linearly with respect to the position within a domain of length L to O(η2). Written in this form it becomes clear that X0 may be interpreted as the macroscopic approximation of X.

Returning to dimensional variables yields the results in the form

(30)

where x=X+χ, X is the macroscopic position, χ=O(L) is the position relative to the macroscopic position, and V and P are the macroscopic particle velocity and acoustic pressure. See Fig. 1 for a summary schematic. These results are valid even for weakly nonlinear waves.

FIG. 1.

Schematic describing the results of the locally valid multiple scales approximation. Consider a pressure wave (the particle velocity may be treated similarly) with characteristic amplitude p0 and minimum time scale T propagating through a one-dimensional heterogeneous domain with a representative volume element length of L. For ηL/cminT1 (cmin is the smallest wave speed present in the domain), the pressure in an arbitrary L-long length of the domain may be written as p(x) with x=X+χ, where X is the center of the L-long length and χ is the position relative to X. The derivation in the main text shows that p(x) is constant with respect to χ to O(ηp0), and is linear with respect to χ to O(η2p0). Thus, the pressure may be written as p(x)=P(X)+χP(X)+O(η2p0)=P(X)+O(ηp0), where P(X) is independent of χ and may be considered to be the macroscopic pressure.

FIG. 1.

Schematic describing the results of the locally valid multiple scales approximation. Consider a pressure wave (the particle velocity may be treated similarly) with characteristic amplitude p0 and minimum time scale T propagating through a one-dimensional heterogeneous domain with a representative volume element length of L. For ηL/cminT1 (cmin is the smallest wave speed present in the domain), the pressure in an arbitrary L-long length of the domain may be written as p(x) with x=X+χ, where X is the center of the L-long length and χ is the position relative to X. The derivation in the main text shows that p(x) is constant with respect to χ to O(ηp0), and is linear with respect to χ to O(η2p0). Thus, the pressure may be written as p(x)=P(X)+χP(X)+O(η2p0)=P(X)+O(ηp0), where P(X) is independent of χ and may be considered to be the macroscopic pressure.

Close modal

It is important to note that the analysis presented above assumed that the inhomogeneities only consist of layers of material with constant mass density and compressibility. Thus, a modified analysis is necessary to account for layers with dynamic material properties, additional physics (e.g., piezoelectricity), or hidden degrees of freedom (e.g., side-wall resonators).

Hamilton's principle is global in nature due in large part to the fact that it is represented as an integral equation, and so naturally the integral itself is important to global homogenization. Consider the spatial integral over some function f(x),

(31)

Any integral may be divided into a summation of sub-integrals that span the same domain as the original integral. Let I be divided into segments of length L, or into M parts such that M=(x1x0)/L:

(32)

See Fig. 2 for a schematic with explanations of the relevant variables. If we identify x0+(m1/2)L with the macroscopic position Xm then the integral may be written as

(33)

where χ=xXm=O(L). The quantity within the square brackets is then identified as the local χ-average of f centered at Xm, and may be written as

(34)

where · is the local χ average, defined as

(35)

As a continuum approximation (which is itself a form of a multiple-scales approximation), the element length L may be approximated as being infinitesimal, and the summation may be approximated by an integral over X:

(36)

This approximation is valid if L is smaller than the smallest characteristic length scale of the macroscopic behavior, i.e., cminT. Thus, the continuum approximation given here is appropriate for η1. Extending this result to larger values of η requires special care (see Ref. 12).

FIG. 2.

Schematic of the globally valid homogenization and continuum approximation. The domain of some multi-scale function f(x) extends from x0 to x1. This domain is partitioned into M subdomains, each L long, with the center of the mth subdomain located at x = Xm. Inside each subdomain there is a local coordinate χ such that in the mth subdomain we may write x=Xm+χ, with L/2<χ<L/2. After the continuum approximation, f(x) is approximated by its local average f, as defined in the text, and the spatial coordinate x is replaced by the macroscopic coordinate X.

FIG. 2.

Schematic of the globally valid homogenization and continuum approximation. The domain of some multi-scale function f(x) extends from x0 to x1. This domain is partitioned into M subdomains, each L long, with the center of the mth subdomain located at x = Xm. Inside each subdomain there is a local coordinate χ such that in the mth subdomain we may write x=Xm+χ, with L/2<χ<L/2. After the continuum approximation, f(x) is approximated by its local average f, as defined in the text, and the spatial coordinate x is replaced by the macroscopic coordinate X.

Close modal

Using the results from Eqs. (30) and (36), Hamilton's principle in Eq. (2) with inhomogeneous material properties may be written approximately as

(37)

where

(38)
(39)
(40)

Since the Lagrange multiplier μ takes whatever form necessary to enforce the constraint, it may be assumed to not be a function of χ and therefore need not be averaged over. The average of the constraint term, however, requires some care. Since the original constraint is just C = 0 for all x, then αC=0 could also have been used, where α is any known function of x. However, α then appears in the constraint average as αC, and is therefore interpreted as a weight. This observation leads to the question, what weight should be chosen? A useful tool to answer this question is the covariance, defined as

(41)

A particularly useful identity of the covariance is that cov(a,b)var(a)var(b), where var(a)=cov(a,a) is the variance of a. Thus, if the variance of a is zero then so is the covariance of a and b, regardless of what b is. With this in mind, consider the constraint of the constitutive equation. The average of this constraint then becomes

(42)

The covariances on the right represent the errors of the continuum approximation, and so it is desirable for them to be as small as possible. From the multiple-scales analysis above we know that the variance of p and v are O(η4) [since p and v are approximated up to O(η2) and var(p) is O(p2)], and so cov(αβ,p)=O(η2) as long as var(αβ)=O(1). However, there is no guarantee that the variance of v/x is small. In fact, throughout the RVE the velocity gradient is proportional to the compressibility times the near-constant pressure. By choosing α = 1 (or any other constant) the final covariance may be set identically to zero, leaving the accuracy of the continuum approximation O(η2). From this example a useful rule is to minimize the variance of any function multiplying the gradient of an acoustic variable.

The relevant variations for Eq. (37) may be written to O(η2) as

(43)
(44)
(45)

Derivatives of the variations of P and V may be moved via integration by parts (boundary terms will be considered below). Thus the integrand of Eq. (37) may be written as

(46)

Eliminating μ from the resulting δV and δP equations yields the homogenized dynamic equation

(47)

Together with Eq. (40) and C=0, the homogenized dynamic equation describes the effective behavior of the system to O(η2), and yields the wave equation

(48)

The effective material properties may then be identified from Eqs. (40) and (47) as

(49)

where ρeff is the effective mass density, βeff is the effective compressibility, and ψeff is a material property that couples the traditional momentum balance equation with the traditional pressure−strain relationship, and has been called the Willis coupling.6,13,14 As has been noted in previous research,15,16 Willis coupling appears only for intrinsically asymmetric microstructures. The effective Willis coupling does not appear in the effective wave equation, as has been noted in previous research,17 and is only manifested at interfaces, making it interesting for applications such as impedance matching and surface waves. It should be noted that Willis coupling is only guaranteed to be absent in one-dimensional effective wave equations; no such guarantee occurs for higher dimensional systems, even if only plane waves are considered.18 

Now consider the boundary terms obtained by integrating Eqs. (43)–(45) by parts. Neglecting temporal boundaries, the boundary terms may be written as

(50)

where U is the macroscopic particle displacement. From Eq. (46) we see that SP+μ/t=O(η), and since χβ=O(η) the coefficients of δP are O(η2) and may be neglected. Therefore, the boundary terms may be written as

(51)

The term μ/t may be substituted to O(η2) with P+(χβ/β)V/t. Following the same logic as found around Eq. (12) we therefore conclude that the macroscopic displacement and the pressure-related quantity

(52)

are continuous at interfaces to O(η2).

It is interesting to note that the macroscopic continuity and dynamic equations may be written in terms of Pˇ, respectively, as

(53)
(54)

These equations lead to a slightly different set of effective material properties, the difference being that the effective Willis coupling ψeff=ρχβ/β is replaced by ψeff=χρ.

The homogenization results derived above are similar to the results of many other methods. However, the present method has the benefit of versatility beyond other methods, as shown in Sec. III.

For the sake of simplicity, all analyses presented below will neglect terms O(η). As a consequence of this simplification, Willis coupling is also neglected.

The homogenization results derived above neglected body forces, but they may be readily accounted for. Consider the case of a one-dimensional heterogeneous acoustic system with weak viscous drag, such that the body force in Eq. (2) may be expressed as19 

(55)

where μvisc is the coefficient of viscosity. Note that this viscous term is actually just the first term of a constitutive equation, and so is only valid up to O(μvisc). The locally valid approximation given above may still be used as long as the viscosity term is negligible over the length scale L, which is the case for η(μvisc/ρcL)max, where the subscript “max” denotes the largest value in the system. In this case, Eq. (37) becomes

(56)

where Ekin,Epot, and C are given as

(57)
(58)
(59)

The homogenization of the left-hand side proceeds as above. The right-hand side is written in terms of gradients of the velocity and the variation of the displacement. Using the local dynamic and constraint equations and integrating-by-parts in time the right-hand side of Eq. (56) may be written in the form

(60)

where no gradients of acoustic variables are present. Then the homogenization procedure given above may be utilized to write

(61)

Thus Eq. (56) becomes

(62)

which leads to the equations

(63)
(64)
(65)

with μ/t=P and U being continuous at interfaces. Combining these three equations leads to the effective dynamic equation

(66)

Since this equation is only valid up to O(μvisc), the quantity 2V/t2 on the right may be substituted with (1/βρ)2V/X2 without changing the order of accuracy, leading to

(67)

where the gradient of the compressibility is assumed to be negligibly small. Thus, the effective homogeneous coefficient of viscosity may be defined as

(68)

The homogenization methods described above are valid for weakly nonlinear systems, such that O(ϵ) terms may be kept without changing the homogenization procedure. This example considers a heterogeneous one-dimensional system without losses that has amplitudes sufficiently high that first-order nonlinear terms must be maintained for accuracy.

The form of the potential energy defined in Eq. (4) is only valid for linear propagation, and is a special case of a more general expression,

(69)

with the pressure and strain being related through Eq. (13), which may be inverted as

(70)

Since the pressure is continuous at interfaces and the strain is not, the most accurate homogenization may be obtained by writing the potential energy entirely in terms of pressure and using Eq. (70) as a constraint equation. In this case and neglecting body forces, Hamilton's principle yields

(71)

Since the strain is isolated from anything except the Lagrange multiplier, this form of Hamilton's principle is appropriate for homogenization, which yields

(72)

where

(73)

is the macroscopic strain. In this case further analysis is performed more naturally in terms of the macroscopic particle displacement U rather than the velocity V, and so we write Hamilton's principle as

(74)

Taking the variation, integrating by parts, collecting like variations and setting the associated coefficients to zero then yields the equations

(75)
(76)
(77)

and μ(1+U/X)δU is continuous at interfaces. These equation may be combined to yield the nonlinear wave equation for the macroscopic displacement:

(78)

The coefficient of the nonlinear term on the right may be identified as the effective coefficient of nonlinearity βnl.20 For homogeneous systems, B may be written in terms of βnl as

(79)

Then, the effective coefficient of nonlinearity may be written as

(80)

It is interesting to compare this result with the effective coefficient of nonlinearity obtained by Everbach et al.,21 which is just βnl,eff=βnlβ2/β2. Note that in the limit of small variance of the compressibility these two predictions coincide. Another note is that while the present analysis is cast in Lagrangian coordinates, the effective material properties do not rely on the coordinate system and are applicable for Eulerian analyses as well.

Incorporating nonlinearity in propagation usually leads to changes in the spectral amplitudes of a signal. In particular, quadratic nonlinearity without sufficient absorption leads to shock formation. Since shocks are characterized by very rapid changes in the acoustic variables, the validity of the multiple scales approximation that is used in the present homogenization scheme comes into question. Signals propagating without microscopic or macroscopic losses will eventually steepen to the point that η/1 and treating the system as an effective medium will become inappropriate. However, for propagation distances much shorter than the macroscopic shock formation distance it may be reasonable to use the effective properties given here.

As a demonstration of the homogenization method, consider a system consisting of 90% water and 10% alcohol placed in periodic layers of length L = 10 cm (see Fig. 3). The constituent materials were chosen due to the ready availability of their traditional and nonlinear material properties. The mass density, compressibility, sound speed, and coefficient of nonlinearity for the water, alcohol, and the effective medium are summarized in Table I. It is interesting to note that the effective mass density and sound speed are very similar to that of water, while the compressibility and, especially, the coefficient of nonlinearity are significantly larger.

FIG. 3.

(Color online) Example of nonlinear homogenization. A periodically heterogeneous system of water (90%) and alcohol (10%) is considered. The length of periodicity is 10 cm, and 100 periods are shown. A high-amplitude Gaussian pulse is initialized at the source at x = 0 on the left and propagates through the system. COMSOL simulations of the actual propagation (Full system) and the propagation assuming effective material properties (Effective system) are presented as a function of distance at several times.

FIG. 3.

(Color online) Example of nonlinear homogenization. A periodically heterogeneous system of water (90%) and alcohol (10%) is considered. The length of periodicity is 10 cm, and 100 periods are shown. A high-amplitude Gaussian pulse is initialized at the source at x = 0 on the left and propagates through the system. COMSOL simulations of the actual propagation (Full system) and the propagation assuming effective material properties (Effective system) are presented as a function of distance at several times.

Close modal
TABLE I.

Relevant material properties for the demonstration of nonlinear homogenization. Properties for water and alcohol obtained from Refs. 22 and 23.

PropertySymbolWaterAlcoholEffective
Mass density (kg/m3ρ 998 790 977 
Sound speed (m/s) c 1480 1150 1421 
Compressibility (1010 Pa−1β 4.6 9.6 5.1 
Coefficient of nonlinearity βnl 3.5 6.3 4.9 
PropertySymbolWaterAlcoholEffective
Mass density (kg/m3ρ 998 790 977 
Sound speed (m/s) c 1480 1150 1421 
Compressibility (1010 Pa−1β 4.6 9.6 5.1 
Coefficient of nonlinearity βnl 3.5 6.3 4.9 

COMSOL Multiphysics 5.4 was used to simulate a signal propagating through the above described system. The source signal is a Gaussian pulse given by

(81)

with p0=29.11 MPa, t0=5L/ceff=0.35 ms, and ceff=1421 m/s is the effective sound speed. The amplitude was chosen such that the effective shock formation distance would be at 8 m. It should be noted that using the Everbach prediction of the effective coefficient of nonlinearity with the above waveform and system would lead to an estimate of the shock formation distance being 8.22 m. The nonlinear propagation was modeled by the second-order Westervelt equation without losses.24 The maximum element size was chosen to be L/40 and the maximum frequency to resolve was chosen to be 10/t0=28.4 kHz.

The COMSOL simulation using the full system and the effective system are shown as a function of distance at several times in Fig. 3. Immediately evident is the nonlinear steepening of the propagating waves, and by 5 ms a shock wave has almost completely formed around the shock formation distance of 8 m. The signals for the full system (solid lines) and the effective system (dashed lines) are very similar for all of the times shown. Of particular note, the shock fronts of the full and effective predictions occur at the same locations for a given time. However, at 4 and 5 ms the full system predictions begin to show small scale high-frequency oscillations that do not appear in the effective system predictions, especially for positions past 5 m. These oscillations are likely due to the fact that shorter wavelengths become important as the waves steepen, making the effective medium approximation become less applicable. This argument may be explored analytically, as shown below.

Since the present analysis is valid to O(η), the quantity η=L/cminT is the important quantity to analyze. The time scale T depends on the time-derivative of the acoustic pressure, which was shown by Muhlestein and Gee to evolve as25 

(82)

By definition p/t=O(p0/T), and so Eq. (82) dictates that

(83)

In the example above, T(0)=t0, and so at the source η=(ceff/cmin)/5=0.25. The approximation is only valid for η<1; this threshold is met at x = 5.17 m, which correlates well with the appearance of the high-frequency oscillations that appear in Fig. 3.

While further quantitative analyses of this example could be given, the qualitative discussion given here is sufficient to demonstrate the efficacy of the present homogenization method.

A one-dimensional acoustic metamaterial homogenization method has been derived and demonstrated. A formal multiple-scales method was used to obtain a locally valid approximations of the acoustic variables. These approximate variables were then substituted into Hamilton's principle, which is global in nature, and then averaged to obtain estimates of the effective material properties of heterogeneous systems. While this homogenization method is approximate in nature, it is also highly versatile. This versatility was demonstrated with two examples: accounting for viscosity, and accounting for finite-amplitude effects.

This method may be considered a framework for additional development. Three possible improvements are mentioned here. First, the present analysis is restricted to one dimension, but there is no fundamental reason that it may not be extended to two and three dimensions. Second, hidden degrees of freedom, such as side-wall resonators for plane wave tubes and mass-in-mass systems, may also be incorporated. Third, additional physics, such as thermal and piezoelectric phenomena, may be accounted for using the method of Lagrange multipliers.

The author would like to thank Carl Hart and Kyle Dunn for excellent discussions of the manuscript at many stages, and the anonymous reviewers for their helpful suggestions and guidance. Permission to publish was granted by the Director, Cold Regions Research and Engineering Laboratory.

1.
J. D.
Eshelby
, “
The determination of the elastic field of an ellipsoidal inclusion, and related problems
,”
Proc. R. Soc. A
241
,
376
396
(
1957
).
2.
A. M.
Baird
,
F. H.
Kerr
, and
D. J.
Townend
, “
Wave propagation in a viscoelastic medium containing fluid-filled microspheres
,”
J. Acoust. Soc. Am.
105
,
1527
1538
(
1999
).
3.
D. V.
Schroeder
,
Introduction to Thermal Physics
(
Addison Wesley Longman
,
Chicago, IL
,
2000
).
4.
R.
Hill
, “
A self-consistent mechanics of composite materials
,”
J. Mech. Phys. Solids
13
,
213
222
(
1965
).
5.
V.
Fokin
,
M.
Ambati
,
C.
Sun
, and
X.
Zhang
, “
Method for retrieving effective properties of locally resonant acoustic metamaterials
,”
Phys. Rev. B
76
,
144302
(
2007
).
6.
M. B.
Muhlestein
,
C. F.
Sieck
,
P. S.
Wilson
, and
M. R.
Haberman
, “
Experimental evidence of Willis coupling in a one-dimensional effective material element
,”
Nat. Commun.
8
,
15625
(
2017
).
7.
P. C.
Andia
,
F.
Costanzo
, and
G. L.
Gray
, “
A Lagrangian-based continuum homogenization approach applicable to molecular dynamics simulations
,”
Intl. J. Sol. Struct.
42
,
6409
6432
(
2005
).
8.
A.
Bedford
,
Hamilton's Principle in Continuum Mechanics
(
Pitman Advanced Publishing Program
,
Boston
,
1985
).
9.
Y.
Fung
and
P.
Tong
,
Classical and Computational Solid Mechanics
(
World Scientific
,
Singapore
,
2001
).
10.
S.
Nemat-Nasser
and
M.
Hori
,
Micromechanics: Overall Properties of Heterogeneous Materials
(
Elsevier
,
Netherlands
,
2013
).
11.
A. H.
Nayfeh
,
Perturbation Methods
(
Wiley
,
New York
,
1973
).
12.
P.
Rosenau
, “
Hamiltonian dynamics of dense chains and lattices: Or how to correct the continuum
,”
Phys. Lett. A
311
,
39
52
(
2003
).
13.
J. R.
Willis
, “
Variational principles for dynamic problems for inhomogeneous elastic media
,”
Wave Motion
3
,
1
11
(
1981
).
14.
G. W.
Milton
and
J. R.
Willis
, “
On modifications of Newton's second law and linear continuum elastodynamics
,”
Proc. R. Soc. A
463
,
855
880
(
2007
).
15.
A. L.
Shuvalov
,
A. A.
Kutsenko
,
A. N.
Norris
, and
O.
Poncelet
, “
Effective Willis constitutive equations for periodically stratified anisotropic elastic media
,”
Proc. R. Soc. A
467
,
1749
1769
(
2011
).
16.
C. F.
Sieck
,
A.
Alù
, and
M. R.
Haberman
, “
Origins of Willis coupling and acoustic bianisotropy in acoustic metamaterials through source-driven homogenization
,”
Phys. Rev. B
96
,
104303
(
2017
).
17.
M.
Muhlestein
and
M. R.
Haberman
, “
Analysis of one-dimensional wave phenomena in Willis materials
,”
Proc. Mtgs. Acoust.
30
,
065017
(
2017
).
18.
M. B.
Muhlestein
,
C. F.
Sieck
,
A.
Alù
, and
M. R.
Haberman
, “
Reciprocity, passivity and causality in Willis materials
,”
Proc. R. Soc. A
472
,
20160604
(
2016
).
19.
A. D.
Pierce
,
Acoustics: An Introduction to Its Physical Principles and Applications
(
Acoustical Society of America
,
Woodbury, NY
,
1989
).
20.
A. N.
Norris
, “
Finite-amplitude waves in solids
,” in
Nonlinear Acoustics
(
Acoustical Society of America
,
Melville, NY
,
2008
), pp.
263
277
.
21.
E. C.
Everbach
,
Z.-m.
Zhu
,
P.
Jiang
,
B. T.
Chu
, and
R. E.
Apfel
, “
A corrected mixture law for B / A
,”
J. Acoust. Soc. Am.
89
,
446
447
(
1991
).
22.
D. T.
Blackstock
,
Fundamentals of Physical Acoustics
(
Wiley
,
Hoboken, NJ
,
2000
).
23.
R. T.
Beyer
, “
The parameter B/A
,” in
Nonlinear Acoustics
(
Acoustical Society of America
,
Melville, NY
,
2008
), pp.
25
39
.
24.
M. F.
Hamilton
and
C. L.
Morfey
, “
Model equations
,” in
Nonlinear Acoustics
(
Acoustical Society of America
,
Melville, NY
,
2008
), pp.
41
63
.
25.
M. B.
Muhlestein
and
K. L.
Gee
, “
Evolution of the temporal slope density function for waves propagating according to the inviscid Burgers equation
,”
J. Acoust. Soc. Am.
139
,
958
967
(
2016
).