## Abstract

An approximate homogenization technique is presented for generally anisotropic elastic metamaterials consisting of an elastic host material containing randomly distributed heterogeneities displaying frequency-dependent material properties. The dynamic response may arise from relaxation processes such as viscoelasticity or from dynamic microstructure. A Green's function approach is used to model elastic inhomogeneities embedded within a uniform elastic matrix as force sources that are excited by a time-varying, spatially uniform displacement field. Assuming dynamic subwavelength inhomogeneities only interact through their volume-averaged fields implies the macroscopic stress and momentum density fields are functions of both the microscopic strain and velocity fields, and may be related to the macroscopic strain and velocity fields through localization tensors. The macroscopic and microscopic fields are combined to yield a homogenization scheme that predicts the local effective stiffness, density and coupling tensors for an effective Willis-type constitutive equation. It is shown that when internal degrees of freedom of the inhomogeneities are present, Willis-type coupling becomes necessary on the macroscale. To demonstrate the utility of the homogenization technique, the effective properties of an isotropic elastic matrix material containing isotropic and anisotropic spherical inhomogeneities, isotropic spheroidal inhomogeneities and isotropic dynamic spherical inhomogeneities are presented and discussed.

## 1. Introduction

The field of acoustic and elastic metamaterials exploits the design of materials with sub-wavelength structure to generate extreme properties on the macroscale. Research on metamaterials can be aided by accurate homogenization models that account for sub-wavelength dynamics. This fact has led to a great deal of research in multiscale effective medium modelling, including elastodynamic homogenization methods. For example, a method that comes from the field of micromechanics is to extend a static homogenization method using Eshelby tensors (see [1] for an introductory discussion) by using dynamic Eshelby tensors [2]. This method has the benefit of being able to account for arbitrarily ellipsoidal inhomogeneities, but does not lend itself easily to the concept of dynamic density or stiffness which are of importance in acoustic and elastic metamaterial research [3]. In another micromechanical approach, Gonnella & Ruzzene [4] were able to approximate the effective in-plane dynamic stiffness and density of two-dimensional periodic lattices using a finite-difference approach with multiple scales. Another approach to dynamic homogenization, single and multiple scattering using multipole expansions of incident and scattered field is commonly used in the field of acoustics and elastodynamics [5]. For example, Baird *et al.* [6] were able to determine analytical expressions for the effective bulk and shear moduli and the effective density of an isotropic elastic sphere with an isotropic elastic coating embedded in an isotropic elastic matrix assuming single scattering. Using a similar approach, Wu *et al.* [7] developed a self-consistent (SC) model to determine the effective dynamic moduli and density for isotropic cylinders embedded in a matrix that shows that negative effective dynamic modulus and density is possible in elastic composites when microscale heterogeneities are resonant. Unfortunately, the mathematical formalisms associated with most scattering-based approaches are restricted to cases of cylindrically or spherically symmetric scatterers and isotropic constitutive materials. For that reason, they are not able to account for anisotropy at either the micro- or the macroscale.

Importantly, most dynamic homogenization methods assume that effective properties take the same form as the properties of the constitutive materials. However, more general constitutive equations are necessary to describe the effective properties of many structures and heterogeneous materials relevant to acoustic and elastic metamaterials (see [8] for a review of some models, and Lakes [9] for a relevant experimental demonstration). Of particular interest for this work, Willis showed that the general macroscopic dynamic response of heterogeneous elastic materials couple stress and momentum due to the existence of ‘polarization’ in both the stress and momentum fields to properly account for the spatial variation of the stiffness and mass density [10–12]. Throughout this work, materials exhibiting this type of coupling are referred to as Willis materials.

The general approach of Willis results in an analytically intractable equation for the effective properties. To address this difficulty, various special cases and approximations have been studied to find explicit expressions for the effective properties. These efforts have primarily been for cases of periodic media, notably those of Nemat-Nasser *et al.* [13–15] and Norris and co-workers [16]. While Sabina & Willis [17] developed an SC homogenization model for random media which formally considers Willis coupling, the final model and examples presented assume symmetry conditions which remove this type of coupling from the resulting effective material properties. Thus, to date there are no models which can estimate macroscopic Willis coupling in media with no long-range order. As some existing and future embodiments of metamaterials involve irregularly placed inhomogeneities, such as so-called ‘soft acoustic metamaterials’ [18], a homogenization technique that permits random distributions of potentially anisotropic inhomogeneities is desired. This paper presents an analytical method that, under certain limiting constraints, allows for homogenization of heterogeneous materials with random distributions of inhomogeneities, including those that display frequency-dependent material properties such as dynamic density, stiffness and Willis coupling, and provides an estimate of the dynamic macroscopic properties.

The structure of this paper is as follows. Section 2 presents a derivation of the homogenization technique, and the associated approximations are explicitly stated. In particular, the important approximation of quasi-static dynamics associated with inhomogeneities embedded in an elastic host material is presented and discussed. Under this approximation, homogenization for arbitrarily anisotropic and dynamic ellipsoidal inhomogeneities may be performed analytically up to a double integral, which may be easily calculated numerically. The resultant effective properties take the form of a coupled constitutive equation, as described by Milton & Willis [12]. Finally, the dilute homogenization technique is extended to higher inhomogeneity volume fractions by employing differential effective medium (DEM) methods. In §3, we demonstrate the utility and limits of the homogenization method with several examples and compare with published methods where applicable. Finally, §4 summarizes the research and draws conclusions about the utility of the model.

## 2. Theoretical development

This section covers the mathematical foundations and approximations of the homogenization technique employed in this work. A brief qualitative overview of the technique is first presented, followed by a rigorous derivation. Whenever an approximation is made in the derivation, a discussion of the physical interpretation and the limits for which the approximation is valid is given.

### (a) Qualitative description

The goal of this homogenization method is to estimate the effective properties of very general elastic metamaterials consisting of inhomogeneities embedded randomly within an elastic matrix material. The resulting model predicts the effective properties of materials with arbitrarily ellipsoidal inhomogeneities that may demonstrate sub-wavelength dynamic behaviour which may yield effective frequency-dependent stiffness, density, and Willis coupling on the microscale (which is the metamaterial approximation). In the present approach, the inhomogeneities may be treated as body force sources within the homogeneous matrix material when subjected to external stimulus. The magnitude of those body force sources is shown to be proportional to the material property contrast between the inhomogeneity and the host material. If the inhomogeneities are very small and well-separated relative to a wavelength, the volume-averaged fields within the inhomogeneities and the volume fraction of the heterogeneities are sufficient to describe their effect on the macroscopic field quantities. Using well-established micromechanical methods [1], the variation of the strain and velocity fields due to the presence of inhomogeneities may be fully described using a Green's function based on the matrix material properties and the geometry of the inhomogeneities. Comparing the resulting strain and velocity fields with a coupled constitutive equation allows effective material properties to be identified. A schematic of the homogenization approach is shown in figure 1. It should be noted that the inhomogeneities may themselves be heterogeneous media in order to generate frequency-dependent effective microscale properties, but these effective properties will not be derived in this paper. While the initial development given here assumes that there is only one kind of inhomogeneity with a single geometric orientation, these results are extended to more higher volume fractions using a DEM theory [19].

While the homogenization technique derived below is only approximate we emphasize that it is very general, and provides the first homogenization method to allow for the estimation of the effective properties of random heterogeneous media displaying macroscopic Willis coupling. Previous homogenization models which account for Willis coupling rely on periodic lattices and have only been demonstrated for one-dimensional systems [13–16], and to the authors' knowledge no numerical homogenization schemes to address random media displaying this type of coupling have been developed. The present method permits the homogenization of fully three-dimensional elastic systems of randomly distributed arbitrarily ellipsoidal inhomogeneities which is valuable for the future design of acoustic and elastic metamaterials.

### (b) Dynamic response of heterogeneous materials

Prior to introducing the homogenization model, it is first necessary to introduce the notation used here to identify tensors for order 0 to 4 as well as the outer product and specific cases of the inner product. Examples of how each order of tensor will be written are shown below.

In the absence of body forces, the dynamic equation of linear elasticity may be written as
** σ** is the Cauchy stress and

_{i}is the divergence over the

*i*th dimension (e.g.

*C*is the stiffness tensor,

*ρ*is the density,

**=**

*ρ**ρ*

**, where**

*δ***is the Kronecker-delta tensor. The stress and momentum depend on the strain and velocity at all times in a material displaying this type of coupling, and this fact is acknowledged by assuming a time harmonic system (**

*δ**e*

^{−iωt}convention) and that the material properties are frequency dependent. Willis coupling may originate from intrinsic asymmetry of the microstructure or from periodicity [15,16], but the latter coupling type is eliminated in the absence of long-range order, as is the case considered here.

Combining the Willis constitutive response with the dynamical equation yields
*M* denotes the matrix material. Clearly, the contrast tensors are zero in the matrix. The spatially varying material properties can then be inserted into equation (2.4) to yield
** G**, is a second rank tensor that satisfies equation (2.8), which is provided in index notation for clarity,

The displacement field may then be implicitly written as
*Ω* which is defined for two arbitrary tensors ** A** and

**as**

*B**T*in equation (2.12) denotes another tensorial transpose which places the first two dimensions at the end of the list of dimensions (e.g.

Now consider a heterogeneous material consisting of two material phases; matrix and inhomogeneity. To mathematically distinguish between the matrix and the inhomogenities, define a Heaviside step function ** X** such that

*I*denotes a property of the inhomogeneity. As the contrast tensors are only non-zero in the inhomogeneities, the convolution volume integrals reduce to integrals over the domain of the inhomogeneities

*Ω*

_{I}. The velocity and displacement gradient fields are then written as

_{2}is the divergence of the second dimension with respect to

*Ω*which is defined for arbitrary

**and**

*A***as**

*B*### (c) Micromechanics

The (exact) velocity and strain fields given in equations (2.14a) and (2.15) are expressed as integral equations whose kernel is an as yet undetermined Green's tensor. In the following, we should note that common micromechanical methods allow these integral equations to be expressed as relatively simple algebraic equations. The micromechanical approximation is that inhomogeneities only interact with other inhomogeneities and affect macroscopic behaviour through volume-averaged quantities. This approximation is valid if the phase of the wave in the matrix varies insignificantly over an inclusion dimension; that is, *k*^{M}*a*≪1, where *k*^{M} is the largest wavenumber in the matrix and *a* is the largest dimension of an inhomogeneity. This condition is also called the quasi-static limit and is of particular interest in the field of metamaterials [20]. Also note that *k*^{eff}*a* must also be small, where *k*^{eff} is the largest wavenumber in the effective material (post-homogenization).

Another consequence of using the quasi-static limit is that the dynamic equation must be approximated for consistency. In performing an order analysis of Willis materials, it is useful to introduce the non-dimensional quantities *k* is the characteristic wavenumber. Thus, ** H**=

**|**

*G**C*|

*a*, we may write equation (2.8) in the dimensionless form

*k*

^{M})

^{2}

*a*

^{2}. Thus, to be consistent with the quasi-static approximation we must neglect the second and third terms, and equation (2.8) reduces to

**is now the elasto-static Green's tensor. As the approximation is only quasi-static, it is equivalent to applying the elastic–viscoelastic correspondence principle and extending it to consider the dispersive effects of subwavelength inclusion dynamics [21]. The quasi-static approximation is of significant importance as it allows one to estimate the effective material properties of interest to acoustic and elastic metamaterials if the frequency-dependent properties of the constituents are known. This is valuable for metamaterial design efforts that may consider heterogeneities with non-spherical geometry and elastic anisotropy of all constituent materials.**

*G*It can be shown that the micromechanical approximation yields
** Y** and

**with the appropriate contraction symbol for °. Equations (2.19) are also known as the mean or average field approximation and are treated in numerous textbooks (e.g. [1,22]). Given this approximation, the volume-averaged strain and velocity fields in a heterogeneous medium where Willis coupling is non-negligible may be written as**

*Z*

*X*^{0}〉 is the homogeneous field

*X*^{0}averaged over the domain of the inclusion and

In order to maintain consistency with the quasi-static approximation an order analysis must be performed on equations (2.20). First, if we define *T*_{1}|∼1/|*C*^{M}|, *T*_{4}|∼*a*^{2}*ω*^{2}/|*C*^{M}|. Then, equations (2.20) may be written in terms of orders of magnitude as
** X**|∼|

*X*^{M}|. As

*k*

^{M}

*aW*

^{M}≪1 and (

*k*

^{M}

*a*)

^{2}≪1 in the quasi-static limit, we find that most of the terms should be neglected at this order of approximation, leaving

*T*

_{1}and

*I*⇔

*I*

_{ijkl}=(

*δ*

_{ik}

*δ*

_{jl}+

*δ*

_{il}

*δ*

_{jk})/2 is the fourth-order symmetric identity tensor.

Because the tensors *A*^{ε}, *A*^{v} map the macroscopic velocity and strain to the volume-averaged microscopic velocity and strain, these four tensors are called ‘localization tensors’ and are a generalized form of the strain localization tensors used in micromechanics. In particular, *A*^{ε} calculates the microscopic strain due to the macroscopic strain, *A*^{v} calculates the microscopic velocity due to the macroscopic velocity. Note that in the quasi-static approximation the macroscopic velocity does not generate a microscopic strain regardless of the coupling tensors. This is due to the fact that the coupling of the velocity to the strain requires a non-negligible phase gradient of the displacement or velocity field on the scale of the inhomogeneity.

In summary, by assuming quasi-static loading in the matrix material (the micromechanical approximation) the exact integral equations for the microscopic displacement and velocity fields may be approximated as algebraic equations which may be solved explicitly in terms of the imposed fields and the localization tensors. While the solution is approximate, it applies to a wide variety of interesting systems, including randomly distributed, fully dynamic inhomogeneities with anisotropic material properties and Willis coupling on the microscale.

### (d) Effective material properties for dilute concentrations

The localization tensors derived in the previous section may be used to predict the effective material properties by extending standard micromechanical homogenization techniques to include the Willis material coefficients. This is done by first assuming
*X*=M, I or eff. Then, using a simple approximation of the overall properties of a two-phase heterogeneous material consisting of inhomogeneities perfectly bonded to the matrix, we may write
*ϕ*≡*Ω*_{I}/[*Ω*_{I}+*Ω*_{M}] is the inhomogeneity volume fraction. The averaged fields 〈** Y**〉

^{M}and 〈

**〉**

*Y*^{I}are the field

**averaged in the matrix and inhomogeneity, respectively, of a representative volume element (RVE). If the inhomogeneities have sufficient spatial separation, interactions between inhomogeneities are negligible, and the RVE may consist of a single inhomogeneity in a matrix. In that case, 〈**

*Y***〉 and**

*ε**Y*〉

^{I}. This is the dilute approximation, which is expressed mathematically as

*k*

^{M}

*a*≪

*k*

^{M}

*L*, where

*a*is the characteristic inhomogeneity size and

*L*is a characteristic distance between inhomogeneities. Since

*ϕ*∼

*a*

^{3}/

*L*

^{3}, the limit where the dilute approximation is strictly valid is

*ϕ*

^{1/3}≪1.

Combining equations (2.28) and (2.27) with equation (2.25), an effective constitutive equation for the heterogeneous material may be found. Under a dilute concentration of heterogeneities in the matrix (*ϕ*≪1), the effective macroscopic field quantities (here denoted with a superscript *g* for ‘global’) may be approximated by the field quantities of a homogeneous matrix; that is, *O*(*ϕ*) [1]. Thus,

It is interesting that in equations (2.34) even if *ρ*^{I}=*ρ*^{I}** δ**,

*ρ*^{M}=

*ρ*

^{M}

**(i.e. the constituent materials satisfy standard constitutive equations), the effective coupling tensors may be non-zero depending on the localization tensors, which take the geometry of the inhomogeneities into account, and associated Green's tensors. Further, equations (2.34) clearly show that coupling may exist even when**

*δ***=0, which is a slight modification from the original observation made by Willis [11], but consistent with later works on the topic [13,16]. It is noteworthy that one observes this result in the quasi-static and dilute limit considered here which implies that Willis coupling may arise from purely local interactions in addition to the existence of non-local interactions associated with an infinite periodic lattice of scatterers, as noted elsewhere [15,16].**

*ρ*### (e) Effective material properties at elevated volume fractions

As stated above, all of these results rely on the assumption of dilute concentrations which is of limited interest. However, well-established methods have been developed to extend this type of development to higher concentrations [1,22]. For example, SC models are based on the idea that if a portion of an effective medium were replaced by the representative volume element of the material the homogenization would cause no change in the effective properties. Another example is DEM theory, which approximates the overall response of non-dilute concentrations of inhomogeneities by iteratively taking a dilute concentration, homogenizing, then treating the resulting effective material as the matrix for another dilute concentration to homogenize. Using these types of methods, the results from §2d may be extended to greater concentrations of inhomogeneities.

Following the derivation shown in [1], the DEM equations which extend the present quasi-static dilute effective material approximation may be written as
*C*^{eff}(*ϕ*=0)=*C*^{M}, *ρ*^{eff}(*ϕ*=0)=*ρ*^{M}, and the localization tensors for the quasi-static case considered here are written as
*O*(*k*^{M}*a*) (or *O*(*k*^{eff}*a*)), and the macroscopic velocity does not generate a microscopic strain. This leads to an estimate of the Willis coupling coefficients where

In addition to higher volume fractions, the present approach can also be extended to multiple inhomogeneity types through iterative methods. The basic idea is that each type or orientation of an inhomogeneity has a certain volume fraction (the sum of which are less than or equal to 1) which may be successively homogenized until all of the full composite is modelled [23]. While this method is straightforward and of interest in the design of metamaterials such as double negative materials, it will not be covered in detail in this paper in the interest of brevity.

## 3. Examples of elastic homogenization

The advantage of the homogenization scheme outlined in §2 is that the localization tensors may be predicted using an analytical model that only requires the evaluation of a simple numerical integral. It is thus an efficient model to explore the effective properties of a heterogeneous material displaying Willis coupling. In this section, we present and briefly discuss four different examples that implement the above homogenization method. In all cases, the heterogeneities are assumed to be ellipsoidal and embedded in a conventional elastic matrix material. The case of isotropic spherical heterogeneities is first presented as a benchmark case, followed by an anisotropic spherical heterogeneity, an isotropic ellipsoidal heterogeneity, and finally an isotropic spherical heterogeneity displaying frequency-dependent mass density and stiffness due to assumed dynamic microstructure. Whenever possible, the predictions from the present model will be compared with other, well-established models, notably, an SC homogenization scheme based on the work of Hill [24], and the Hashin–Shtrikman (HS) bounds for isotropic, spherical inclusions. As there are no published models, analytical or numerical, which account for Willis coupling in random media, the predictions for the coupling tensors will be presented alone.

### (a) Isotropic spherical inhomogeneities

Because there are numerous homogenization techniques which can account for isotropic spherical inhomogeneities in an isotropic matrix, this case may be treated as a benchmark for the present homogenization method. We consider a matrix with shear modulus *μ*^{M}=1 GPa, bulk modulus *K*^{M}=3 GPa, and density *ρ*^{M}=700 kg m^{−3} containing 0.1 mm radius spherical inhomogeneities with shear modulus *μ*^{I}=10 GPa, bulk modulus *K*^{I}=30 GPa and density *ρ*^{I}=2000 kg m^{−3} (heavy and stiff inhomogeneities). Neither the matrix nor the inhomogeneity has any Willis coupling, because the only isotropic third-order tensor is the null third-order tensor. The frequency is assumed to be 50 Hz.

Figure 2 shows the effective plane wave modulus as a function of inhomogeneity volume fraction using the W-DEM approach, equations (2.35) with equations (2.36), with the associated HS bounds and an SC prediction [25,26] based on the work of Hill [24], shown for comparison. The relative error of the W-DEM prediction relative to the SC prediction is 0.3% at 10% volume fraction and is 1.7% at 20% volume fraction. As SC and W-DEM are different homogenization models, it is not expected that they would predict the same effective properties. What is expected is that the two predictions should be similar and to fall within the HS bounds, as is seen in this case.

Because both matrix and inhomogeneity materials are isotropic and the geometry is also isotropic, we expect the effective material properties to also be isotropic and therefore uncoupled. Indeed, this is the case, as seen by the effective material properties summarized in table 1, which gives the matrix, inhomogeneity, W-DEM and SC predicted effective properties, along with the HS bounds for comparison at 5% volume fraction.

### (b) Anisotropic spherical inhomogeneities

Consider a spherical inhomogeneity with anisotropic material properties. For convenience, only transverse isotropy (symmetric about the 1 axis), which allows for 10 independent material constants (5 stiffnesses+2 densities+3 couplings), will be considered. The matrix and inhomogeneity properties are summarized in table 2, along with the effective properties predicted by the W-DEM. The matrix properties are the same as those presented in table 1, the frequency is 50 Hz, the inhomogeneity radius is 0.1 mm, and the inhomogeneity volume fraction is 5%. Figure 3 shows the W-DEM and SC predictions for the stiffness moduli *C*_{1111} and *C*_{2222} (*a*) and the W-DEM prediction for the independent coupling moduli *S*_{111}, *S*_{221} and *S*_{122}.

As expected, the effective material properties are also transversely isotropic and symmetric about the *x*_{1}-axis. All of the predicted material stiffnesses are closer to the Reuss average (lower bound) than the Voigt average (upper bound), which is similar to the isotropic case considered in §3a using the HS bounds. The W-DEM and SC predictions for the stiffness moduli are in table 2, but follow different paths in figure 3, which means that the W-DEM and SC models are nearly the same for small volume fractions and close to one, but diverge somewhat at intermediate volume fractions. The W-DEM prediction of the coupling moduli are much smaller than the Voigt averages, following the pattern of the predicted stiffness moduli. The similarity of the evolution of the coupling moduli to the stiffness moduli is seen at all volume fractions in figure 3.

### (c) Isotropic ellipsoidal inhomogeneities

In addition to material anisotropy, the present homogenization technique can account for anisotropy resulting from the presence of microscale ellipsoidal inhomogeneities. For convenience only spheroids, which are ellipsoids with two equal radii, are considered. This form of inhomogeneity will lead to transversely isotropic effective material properties. In particular, the cases of prolate spheroidal (needle-like) inhomogeneities and oblate spheroidal (penny-like) inhomogeneities will be considered. For both cases, the matrix and inhomogeneity material properties are assumed to be isotropic and the same as those in §3a, and again the frequency is 50 Hz.

The predicted material properties for these two cases are given in table 3. The numerical values are rounded to two digits, as the numerical error associated with the W-DEM prediction precluded any greater accuracy given the symmetry requirements on the material properties. Using an oblate spheroidal inhomogeneity causes the material to be more compliant in the 1 direction than using a prolate spheroidal inhomogeneity, though the stiffness of the inhomogeneity material causes both to be stiffer than the matrix material. Conversely, the stiffness in the 2 and 3 directions is greater for the oblate case than for the prolate case. Interestingly, no coupling was generated in either the prolate or oblate cases because the inhomogeneities are symmetric about the principal axes, and the density in both cases is simply the Voigt average.

The W-DEM predictions for the spheroidal inhomogeneities are compared with the SC predictions in figure 4. Again, the two predictions follow very similar trends for all volume fractions.

### (d) Dynamic mass and stiffness

In addition to dealing with material or geometric anisotropy, the W-DEM approach can account for dynamic inhomogeneities by homogenizing frequency by frequency as long as the long-wavelength limit holds in the effective and constitutive media. For example, consider the isotropic spherical inhomogeneities in §3a with 20% volume fraction where the material properties of the heterogeneities are frequency dependent. This can be simply modelled by multiplying the property of interest by a frequency-dependent factor *F*(*ω*) representing dispersion. For illustrative purposes, consider the case where *F*(*ω*) is chosen to be a Lorentzian distribution defined as
*ω*_{0} is some (very small) characteristic frequency and *Γ* is a damping factor. The function *F*(*ω*) has the properties *F*(0)=1 and *Γ*=0, *F*(*ω*_{0})=0 (see figure 5 for a plot of *F*(*ω*) with *Γ*=10^{−2}, 10^{−1} and 1). The maximum of the real part of *F*(*ω*) occurs at *F*(*ω*) follows the form given by Fang *et al.* [27] which represents the effective stiffness in a one-dimensional periodic array of Helmholtz resonators.

For this example, the stiffness and mass density tensors of the inhomogeneity are both multiplied by *F*(*ω*). Varying *C*^{I} in this manner is equivalent to holding Poisson's ratio constant and varying the shear or bulk modulus. Evaluation of the W-DEM with uncoupled constitutive materials yields an effective stiffness tensor that is insensitive to the constitutive mass densities, and the effective mass density tensor is insensitive to the constitutive stiffness tensors. Therefore, both the effects of frequency-dependent stiffness and mass density may be analysed at the same time, which is done here.

The W-DEM prediction of the dispersion curves for *μ*^{eff} and *ρ*^{eff} is shown in figure 6. It is interesting to note the strong dispersion of *μ*^{eff} occurs at *ω*/*ω*_{0}≈1, while the frequency range of strong dispersion of *ρ*^{eff} occurs at *ω*/*ω*_{0}≈1.7. This is due to the natures of the localization tensors *A*^{ε} and *A*^{v} in the W-DEM. Since *A*^{v}≈** δ** (equation (2.36d)) and

*F*(

*ω*) (which is at

*ω*=1.73

*ω*

_{0}in the undamped case). However, A

^{ε}≈[I−T:ΔC]

^{−1}(equation (2.36a)), so the strongly dispersive response of the effective stiffness given in equation (2.34a) will occur when

*T*:Δ

*C*≈

*I*. It should be noted that since the resonance condition depends upon

*C*

^{M}in the dilute approximation, the implicit nature of the W-DEM method makes explicit analysis of this phenomenon difficult for significant volume fractions.

## 4. Summary

A homogenization technique for arbitrarily anisotropic elastic heterogeneous materials with no long-range order that exhibit Willis coupling has been derived and demonstrated. Explicit prediction of the effective material properties required using a micromechanical approximation (*k*^{M}*a*≪1, where *a* is the characteristic dimension of the inhomogeneity) to turn the integral equations of Willis into algebraic equations that are presented for a biphase heterogeneous material with non-interacting inhomogeneities (inhomogeneity volume fraction *ϕ*≪1). Using a standard micromechanical homogenization approach to approximate the macroscopic stress, strain, momentum and velocity fields, the effective material properties may be approximated using localization tensors. As the approximation neglects effects *O*(*k*^{M}*a*), only the strain localization is accurately predicted, and so only one of the Willis coupling tensors, *k*^{M}*a* is left as future work.

The homogenization technique was used to predict the effective material properties of various systems using a DEM approach to extend the results to all 0≤*ϕ*≤1. The results are very similar in magnitude and volume fraction dependence to an established SC homogenization scheme and falls within the well-established HS bounds when Willis coupling is neglected. It was also found in the cases shown here that the presence of Willis coupling in the inhomogeneities does not significantly affect the effective stiffness. It was demonstrated that Willis coupling is not the result of anisotropy emerging from the presence of non-spherical inhomogeneities with three axes of symmetry, as isotropic ellipsoidal inhomogeneities embedded in an isotropic matrix do not generate any coupling. However, the presence of Willis coupling in an inhomogeneity does generate Willis coupling in the effective material properties. These findings are consistent with recent work suggesting that Willis coupling is the result of long range order and a lack of ‘reflective symmetry’ on the microscale [15,16]. Finally, the ability to calculate effective dynamic stiffness and density when dynamic inhomogeneities having known, frequency-dependent, material properties are present has been shown. The behaviour of the effective dynamic stiffness is shown to behave differently than the effective density, even when the inhomogeneity stiffness and density have similar microscale dynamics. This homogenization technique may be used to design acoustic and elastic metamaterials with random microstructure in cases where the matrix material and/or heterogeneities display frequency-dependent stiffness, mass density and Willis coupling as a result of microscopic inclusion structure. This model may be expanded to investigate the influence of multiple types of heterogeneities as well as the effects of orientation distribution in a random heterogeneous medium.

## Authors' contributions

M.R.H. conceived the original homogenization method presented in the paper and provided valuable guidance. M.B.M. developed, expanded and implemented the homogenization method. Both authors gave final approval for publication.

## Competing interests

We have no competing interests.

## Funding

This work was supported by ONR through MURI grant no. N00014-13-1-0631.

## Acknowledgements

The work for this paper was performed independently by the authors.

## Appendix A. Evaluation of *T*_{1} and T 3 tensors for ellipsoidal inhomogeneities

Throughout these derivations, extensive use of index notation will be made for the sake of clarity. Also, all material properties will be assumed to be matrix properties unless otherwise specified. In preparation to use spectral methods to analyse the tensors *T*_{1} and *k*_{i}=*kχ*_{i} is the wavevector, *k* is the wavenumber, and *χ*_{i} is the unit direction-cosine vector. Note the transforms have the property

**(a) Derivation of **

The Fourier transform of the quasi-static Green's function equation given in equation (2.18) yields *P*_{ik} = *χ*_{l}*χ*_{j}*C*_{ijkl}.

The tensor *T*_{1} as defined in equation (2.21a) may be expanded as *Ω*_{I} is the inhomogeneity volume. The integrals over physical space may be performed exactly, as can the resulting integral over *k*. However, the remaining integrals over *θ* and *ϕ* cannot be readily integrated analytically. If we define
*T*_{1})_{ijkl} may be written *P*_{ij}|∼|*C*_{ijkl}|, |*T*_{ijkl}|∼|*C*_{ijkl}|^{−1} and is independent of inhomogeneity scale.

**(b) Derivation of **

The tensor *G*_{ik,j} for its Fourier transform and performing similar integrals as in the analysis of

**(c) Extension to ellipsoidal inhomogeneities**

Consider an ellipsoidal inhomogeneity defined such that the outer boundary of the inhomogeneity is expressed as
*X*_{i} such that *x*_{i}=*Φ*_{ij}*X*_{j}, where
*K*_{i} should be defined such that *a*_{1} and *x*_{i}*k*_{i}=*X*_{i}*K*_{i}. The Jacobian determinant of the spatial transformation is *Ω*_{k} is the *k*-space volume. Recall *χ*_{i} is the vector of direction cosines of *k*_{i}. Thus, *χ*_{i}′ is the vector of direction cosines of *K*_{i}. Defining *P*′_{ik}≡*χ*′_{l}*χ*′_{j}*C*′_{ijkl}, where

**(d) Analytical forms of the tensors for an isotropic matrix with spherical inhomogeneities**

An important case to consider is an isotropic matrix with spherical inhomogeneities. In this simple case, analytical solutions may be found for *C*_{ijkl}=*λδ*_{ij}*δ*_{kl}+*μ*(*δ*_{ik}*δ*_{jl}+*δ*_{il}*δ*_{jk}). For this case, the tensor *P*_{ik} may be written as

- Received June 3, 2016.
- Accepted August 2, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.