Search tips
Search criteria 


Logo of procaThe Royal Society PublishingProceedings AAboutBrowse by SubjectAlertsFree Trial
Proc Math Phys Eng Sci. 2016 August; 472(2192): 20160438.
PMCID: PMC5014117

A micromechanical approach for homogenization of elastic metamaterials with dynamic microstructure


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.

Keywords: metamaterials, micromechanics, elastic homogenization, ellipsoidal inhomogeneity

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 [1012]. 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. [1315] 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].

Figure 1.
A schematic of the homogenization method. The scattered field (red dashed lines) from the inhomogeneities (yellow ellipses) due to an incident field (red solid lines) may be considered point sources in a matrix (grey background). Effective material properties ...

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 [1316], 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.

ϕ(unadorned)zeroth order (scalar),u(arrow)first order (vector),σ(bold)second order,S(calligraphic)third order,C(sans-serif)fourth order.

For simplicity, the coordinate system is assumed to be Cartesian and standard index notation will be used when convenient. There are various types of inner product, or contractions. The dot notation will be used, such that a single dot sums over the product of the last index of the left tensor, and the first index of the right tensor (e.g.  ABAijBjkl). Similarly, two dots sum over the last two indexes of the left tensor and the first two indexes of the right tensor (e.g.  A:BAijBijk). We are now in a position to introduce the homogenization model.

In the absence of body forces, the dynamic equation of linear elasticity may be written as


where σ is the Cauchy stress and p is the momentum density. The operator divi is the divergence over the ith dimension (e.g. div2σσij,j), and over-set dots denote partial differentiation with respect to time. These equations describe the material response to an imposed deformation or traction and must be supplemented with constitutive equations. The standard choices for linear elasticity are


where C is the stiffness tensor, ε=(u+(u)T)/2 is the linear strain tensor, ρ is the density, u is the displacement field, and u is the displacement gradient field (i.e. uui,j). The raised ‘T’ denotes a transpose which switches the order of the first two dimensions of the tensor (e.g. Gij,kT=Gji,k). However, in the more general case of the effective dynamic response of heterogeneous materials, a more general set of constitutive equations known as the Willis equations [10] should be employed. The constitutive relations for the stress and momentum fields in a Willis material are


where S and S~ are effective material properties that couple the macroscopic stress and velocity fields and the macroscopic momentum and strain fields. This response is the result of microscale inertia that may result from macroscopically imposed strain when dynamic microstructure is present [12]. Note that the Willis equations reduce to the standard constitutive equations for S=S~=0 and ρ=ρδ, 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


An elegant way to deal with the spatial dependence of the material properties in piece-wise heterogeneous media is to introduce contrast tensors defined as


for X{C,S,S~,ρ}, where the superscript 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







As the notation implies, the material contrast terms may be treated as force source distributions in a homogeneous matrix material. The effect of these sources on the overall response of the material may therefore be described using the Green's function approach described later. Because the displacement field is a vector quantity and the source terms are vector quantities, the Green's function, G, is a second rank tensor that satisfies equation (2.8), which is provided in index notation for clarity,


The above equation is for point-source-driven wave motion in a general Willis material, where  G(x) is the frequency domain Green's function and δ(x) is the three-dimensional Dirac delta function. The Green's function gives the displacement field at the position x due to a point force source at the origin (or at a source location x if the arguments of  G(x) and δ(x) are replaced by xx). Furthermore,  G(x) is material, geometry and frequency dependent, in general.

The displacement field may then be implicitly written as


where u0(x) is the solution to the homogeneous boundary-value problem. The binary operator Ω is a tensorial spatial convolution over the domain Ω which is defined for two arbitrary tensors A and B as


Note that for time harmonic systems


and so


Thus, the integral equation for the velocity may be written as


Another property of convolutions, ( AΩ B)=( A)Ω B= AΩ( B), then allows us to write


The raised sanserif 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. Gij,kT=Gki,j).

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) which is equal to zero if x is in the matrix and equal to one if x is in the inhomogeneity. Now define ΔX such that δ X(x)=θ(x)Δ X=θ(x)( XI XM) for  X{C,S,S~,ρ}, where the superscript 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




where div′2 is the divergence of the second dimension with respect to x. Note that the solution to the homogeneous boundary-value problem u0(x) is still accounted for when xΩI. By an extension of the properties of derivatives of convolutions to tensorial spatial convolutions, the velocity and displacement gradient fields become




where the binary operator Ω is another tensorial spatial convolution over the domain Ω which is defined for arbitrary A and B as


Using equation (2.14b), the strain field can then be written as


This formulation of the velocity and strain is analogous to the formulation given by Sabina & Willis [17].

(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, kMa[double less-than sign]1, where kM 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 keffa must also be small, where keff 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 W=S/Z0 and W~=S~/Z0, where Z0=|C||ρ| is the characteristic specific acoustic impedance of the material and |[center dot]| denotes the characteristic amplitude of the argument. The interpretation of W and W~ becomes apparent when considering the order of magnitude of the terms in the constitutive equations (2.3):


where k is the characteristic wavenumber. Thus, |W| and |W~| represent the relative importance of Willis coupling in the constitutive equations. By reciprocity and passivity, we also know that |W~|=|W| [15,16]. If we define y=x/a and H=G|C|a, we may write equation (2.8) in the dimensionless form


where partial differentiation is now with respect to y and we have assumed kω|C|/|ρ|, which is equivalent to assuming |W|21. Note that in the case of strong coupling, or |W|21, the wavenumber is approximately k=i|W|ω|C|/|ρ|. The characteristic amplitudes of the three coefficients in equation (2.17) are, in order from left to right, 1, kMa|W|, and (kM)2a2. Thus, to be consistent with the quasi-static approximation we must neglect the second and third terms, and equation (2.8) reduces to


where G 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.

It can be shown that the micromechanical approximation yields


for arbitrary tensors Y and Z 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




where left angle bracketX0right angle bracket is the homogeneous field X0 averaged over the domain of the inclusion and









The expressions for N(x) and M(x) have been provided in index notation in equations (2.22) for clarity.

In order to maintain consistency with the quasi-static approximation an order analysis must be performed on equations (2.20). First, if we define h=u/a, then we may write xu=yh and u˙=iωah. Also, note that |T1|~1/|CM|, |T2||T3|aω/|CM|, and |T4|~a2ω2/|CM|. Then, equations (2.20) may be written in terms of orders of magnitude as




where we have assumed |ΔX|~|XM|. As kMaWM[double less-than sign]1 and (kMa)2[double less-than sign]1 in the quasi-static limit, we find that most of the terms should be neglected at this order of approximation, leaving




Expressions for T1 and T3 may be written down analytically for isotropic spherical heterogeneities and with some numerical integration for ellipsoids in a quasi-static system (see appendix A). Equations (2.24) can then be inverted to obtain an expression for the volume-averaged strain and velocity fields in the heterogeneity as a function of the macroscopic field values









where I[left and right double arrow ]Iijkl=(δikδjl+δilδjk)/2 is the fourth-order symmetric identity tensor.

Because the tensors Aε, Aεv, Avε and Av 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 strain due to the macroscopic velocity, Avε calculates the microscopic velocity due to the macroscopic strain and Av 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




where 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


for  Y{σ,p,ε,v}, where the superscript eff denotes the effective, or macroscopic, property and ϕ[equivalent]ΩI/[ΩI+ΩM] is the inhomogeneity volume fraction. The averaged fields left angle bracketYright angle bracketM and left angle bracketYright angle bracketI are the field Y 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, left angle bracketεright angle bracket and u˙ in equation (2.25) may be used to represent left angle bracketY right angle bracketI. This is the dilute approximation, which is expressed mathematically as kMa[double less-than sign]kML, where a is the characteristic inhomogeneity size and L is a characteristic distance between inhomogeneities. Since ϕ~a3/L3, the limit where the dilute approximation is strictly valid is ϕ1/3[double less-than sign]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 (ϕ[double less-than sign]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, ugu0. This is equivalent to the classical Eshelby method and is valid up to O(ϕ) [1]. Thus,


and also


where we have substituted the matrix strain and velocity similar to the substitution of the matrix stress in equation (2.29), which combine to give


Substituting the field quantities within the inhomogeneities for the general quantities using the localization tensors in equations (2.25), one can write


Similarly, the macroscopic momentum may be written as


Note that the macroscopic stress and momentum fields are related to both the macroscopic strain and velocity fields via the localization tensors, which are linear operators. Thus, by matching equations (2.32) and (2.33) with the assumed form of the effective medium provided in equation (2.27), the effective material satisfies a set of effective Willis equations where





As should be expected, the effective material properties are tacitly frequency dependent based on the dispersion relations of the constitutive materials but independent of the imposed fields.

It is interesting that in equations (2.34) even if SI=S~I=SM=S~M=0 and ρ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 SM=SI=S~M=S~I=0 and Δρ=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





with the conditions Ceff(ϕ=0)=CM, Seff(ϕ=0)=SM, S~eff(ϕ=0)=S~M and ρeff(ϕ=0)=ρM, and the localization tensors for the quasi-static case considered here are written as





These equations may be solved using numerical methods. The results shown in §3 are found using this DEM theory with a fourth-order Runge–Kutta method. As discussed above, the fact that ADEMεv=0 comes from neglecting terms O(kMa) (or O(keffa)), and the macroscopic velocity does not generate a microscopic strain. This leads to an estimate of the Willis coupling coefficients where (S~eff)TSeff. This inequality, which violates reciprocity [15,16], is a consequence of the quasi-static approximation. Despite this fact, the prediction of S~eff provided by this model captures the predominant microscale physics that lead to the macroscale properties of interest and is therefore used as the estimate of macroscopic coupling. This is the case because the model accurately localizes strain, but not velocity, and assumes that dynamic density, stiffness and Willis coupling exists on the microscale (matrix and heterogeneity). Therefore, the macroscale strain and velocity will localize to generate the approximate microscale momentum that depends on both the microscale strain and velocity. Specifically, equations (2.25) and the constitutive relations yield the following expression for the momentum in the heterogeneities:


This momentum density is clearly dependent on both the macroscopic strain and velocity, despite the quasi-static approximation. This localized momentum will, in turn, give rise to macroscopic momentum that depends on the volume fraction, geometry and anisotropy of the coupled inclusions along with the material properties of the matrix through the evaluation of the relations in equation (2.35). It therefore provides a good approximation of S~eff, which is used as the effective Willis coupling coefficient from this point forward. Despite these obvious simplifications, it is important to note that the model provides an estimate of the effective Willis coupling for a class of materials that has not yet been accurately modelled: random anisotropic heterogeneous media demonstrating Willis coupling.

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 KM=3 GPa, and density ρM=700 kg m−3 containing 0.1 mm radius spherical inhomogeneities with shear modulus μI=10 GPa, bulk modulus KI=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.

Figure 2.
Predicted effective plane wave modulus (C1111) for an isotropic matrix with spherical isotropic inhomogeneities using the Willis differential effective medium approach (W-DEM). The HS bounds and the predicted plane wave modulus using an SC approach are ...

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.

Table 1.
List of material properties for an isotropic matrix containing isotropic spherical inhomogeneities with a 5% volume fraction. The predicted material properties using a self-consistent (SC) prediction and the upper and lower HS bounds are shown for comparison. ...

(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 C1111 and C2222 (a) and the W-DEM prediction for the independent coupling moduli S111, S221 and S122.

Figure 3.
(a) Effective plane wave moduli (C1111 and C2222) for an isotropic matrix with identically oriented spherical transversely isotropic inhomogeneities using the W-DEM approach and the generalized SC methods. (b) Effective coupling modulus S111 of the same ...
Table 2.
List of material properties for an isotropic matrix containing spherical inhomogeneities with transversely isotropic material properties. The volume fraction of heterogeneity is 5%. The material properties predicted using the Reuss and Voigt averages ...

As expected, the effective material properties are also transversely isotropic and symmetric about the x1-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.

Table 3.
List of material properties for an isotropic matrix with isotropic prolate or isotropic oblate spheroidal inhomogeneities. The axis of symmetry is the x1, and the ratio of the long to short radius is 10. A volume fraction of 5% is assumed.

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.

Figure 4.
Effective plane wave moduli (C1111 and C2222) for an isotropic matrix with (a) prolate spheroidal (needle-like) or (b) oblate spheroidal (penny-like) isotropic inhomogeneities using the W-DEM approach and the generalized SC methods. (Online version in ...

(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


where ω0 is some (very small) characteristic frequency and Γ is a damping factor. The function F(ω) has the properties F(0)=1 and F()=3, and for Γ=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 ω=ω03+Γ3. The choice of F(ω) follows the form given by Fang et al. [27] which represents the effective stiffness in a one-dimensional periodic array of Helmholtz resonators.

Figure 5.
(a,b) Plot of the characteristic dispersion relation F(ω) as a function of normalized frequency.

For this example, the stiffness and mass density tensors of the inhomogeneity are both multiplied by F(ω). Varying CI 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 Av in the W-DEM. Since Avδ (equation (2.36d)) and ΔS~=0, the effective density given in equation (2.34d) is approximately just a mixture law, and so the strong dispersion of the density will occur at the resonance of 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 TCI. It should be noted that since the resonance condition depends upon CM in the dilute approximation, the implicit nature of the W-DEM method makes explicit analysis of this phenomenon difficult for significant volume fractions.

Figure 6.
Real (a,c) and imaginary (b,d) parts of the W-DEM prediction of the effective shear modulus (μeff, (a,b)) and mass density (ρeff, (c,d)) for an isotropic matrix with dynamic isotropic spherical inhomogeneities as a function of frequency ...

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 (kMa[double less-than sign]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 ϕ[double less-than sign]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(kMa), only the strain localization is accurately predicted, and so only one of the Willis coupling tensors, S~eff, is likely to be well approximated by the homogenization scheme. Extension of the present method to higher orders of kMa 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.


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

Appendix A. Evaluation of T1 and T3 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 T1 and T3, define the three-dimensional spatial Fourier transform (at a single frequency) with the pair

A 1

where ki=i is the wavevector, k is the wavenumber, and χi is the unit direction-cosine vector. Note the transforms have the property G~ij,l=iklG~ij=ikχlG~ij. For later use, note

A 2

(a) Derivation of T1 for spherical inhomogeneities

The Fourier transform of the quasi-static Green's function equation given in equation (2.18) yields k2χlχjCijkl0G~km+δim=0, which may be solved for Green's function as G~ij=k2Pij1, where Pik = χlχjCijkl.

The tensor T1 as defined in equation (2.21a) may be expanded as (T1)ijkl=Tijkl(1)+Tjikl(1), where

A 3

Ω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

A 4

which must be integrated numerically, the tensor Tijkl(1) may be written as

A 5

Therefore, the full tensor (T1)ijkl may be written Tijkl=(Kijkl(4)+Kjikl(4))/8π. This form is equivalent to previously published forms [28,29]. As |Pij|~|Cijkl|, |Tijkl|~|Cijkl|−1 and is independent of inhomogeneity scale.

(b) Derivation of T3 for spherical inhomogeneities

The tensor T3 as defined in equation (2.21b) may be written in index notation as

A 6

Substituting Gik,j for its Fourier transform and performing similar integrals as in the analysis of Tijkl(1), we may write (T3)ijk=3aωKijk(3)/16π2, where

A 7

Thus, |T3|aω/|Cijkl|. Note that since Gij,k(xx)=Gij,k(xx), the value of T3 must be zero for centre-symmetric geometries, and so Kijk(3)=0.

(c) Extension to ellipsoidal inhomogeneities

Consider an ellipsoidal inhomogeneity defined such that the outer boundary of the inhomogeneity is expressed as

A 8

Define a new coordinate system Xi such that xi=ΦijXj, where

A 9

This suggests Ki should be defined such that ki=Φij1Kj. In this coordinate system, the inhomogeneity is a sphere of radius a1 and xiki=XiKi. The Jacobian determinant of the spatial transformation is |Φij|=a2a3/a12, and the Jacobian determinant of the wavevector transformation is |Φij1|=a12/a2a3. Then Tijkl(1) may be written as

A 10

where Ωk is the k-space volume. Recall χi is the vector of direction cosines of ki. Thus, χi=k1ki=k1Φij1Kj=(K/k)Φij1χj, where χi′ is the vector of direction cosines of Ki. Defining Pik[equivalent]χlχjCijkl, where CijklΦmi1Φnj1Φpk1Φql1Cmnpq, such that Pik=(k/K)2Φpi1Φqk1Ppq, results in

A 11

which becomes Tijkl(1)=Φiα1Φjβ1Φkγ1Φlδ1Kαβγδ(4)/8π, where

A 12

Following similar steps, (T3)ijk=3a1ωΦiα1Φjβ1Φkγ1Kαβγ(3)/16π2, where

A 13

(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 Kijkl(4) and Kijk(3). The stiffness tensor for an isotropic material may be written as Cijkl=λδijδkl+μ(δikδjl+δilδjk). For this case, the tensor Pik may be written as

A 14

and the inverse may be written as Pij1=δij/μχiχj/2μ(1ν). Note Pji1=Pij1 is a symmetric tensor. The needed integrals are then written

A 15


A 16

Then for an isotropic matrix with spherical inhomogeneities T3=0 and

A 17

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.


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


1. Qu J, Cherkaoui M 2006. Fundamentals of micromechanics of solids. New York, NY: John Wiley & Sons Inc.
2. Michelitsch TM, Gao H, Levin VM 2003. Dynamic Eshelby tensor and potentials for ellipsoidal inclusions. Proc. R. Soc. Lond. A 459, 863–890. (doi:10.1098/rspa.2002.1054)
3. Cummer SA, Christensen J, Alù A 2016. Controlling sound with acoustic metamaterials. Nat. Rev. Mat. 1, 16001 (doi:10.1038/natrevmats.2016.1)
4. Gonella S, Ruzzene M 2008. Homogenization and equivalent in-plane properties of two-dimensional periodic lattices. Int. J. Solids Struct. 45, 2897–2915. (doi:10.1016/j.ijengsci.2014.02.022)
5. Waterman PC, Truell R 1961. Multiple scattering of waves. J. Math. Phys. 2, 513–537. (doi:10.1063/1.1703737)
6. Baird AM, Kerr FH, Townend DJ 1999. Wave propagation in a viscoelastic medium containing fluid-filled microspheres. J. Acoust. Soc. Am. 105, 1527–1538. (doi:10.1121/1.426692)
7. Wu Y, Lai Y, Zhang Z 2007. Effective medium theory for elastic metamaterials in two dimensions. Phys. Rev. B. 76, 205313 (doi:10.1103/PhysRevB.76.205313)
8. Del Vescovo D, Giorgio I 2014. Dynamic problems for metamaterials: review of existing models and ideas for further research. Int. J. Eng. Sci. 80, 153–172. (doi:10.1016/j.ijengsci.2014.02.022)
9. Lakes R. 1995. On the torsional properties of single osteons. J. Biomech. 28, 1409–1410. (doi:10.1016/0021-9290(94)90260-7) [PubMed]
10. Willis JR. 1981. Variational principles for dynamic problems for inhomogeneous elastic media. Wave Motion 3, 1–11. (doi:10.1016/0165-2125(81)90008-1)
11. Willis JR. 1997. Dynamics of composites. In Continuum micromechanics (ed. P Suquet), CISM Lecture Notes, vol. 495, pp. 265–290. Berlin, Germany: Springer.
12. Milton GW, Willis JR 2007. On modifications of Newton's second law and linear continuum elastodynamics. Proc. R. Soc. A 463, 855–880. (doi:10.1098/rspa.2006.1795)
13. Nemat-Nasser S, Srivastava A 2011. Overall dynamic constitutive relations of layered elastic composites. J. Mech. Phys. Solids 59, 1953–1965. (doi:10.1016/j.jmps.2011.07.008)
14. Srivastava A, Nemat-Nasser S 2011. Universal theorems for total energy of the dynamics of linearly elastic heterogeneous solids. Mech. Mater. 43, 913–917. (doi:10.1016/j.mechmat.2011.06.011)
15. Srivastava A, Nemat-Nasser S 2012. Overall dynamic properties of three-dimensional periodic elastic composites. Proc. R. Soc. A 468, 269–287. (doi:10.1098/rspa.2011.0440)
16. Norris AN, Shuvalov AL, Kutsenko AA 2012. Analytical formulation of three-dimensional dynamic homogenization for periodic elastic systems. Proc. R. Soc. A 468, 1629–1651. (doi:10.1098/rspa.2011.0698)
17. Sabina FJ, Willis JR 1988. A simple self-consistent analysis of wave propagation in particulate composites. Wave Motion 10, 127–142. (doi:10.1016/0165-2125(88)90038-8)
18. Forrester DM, Pinfield VJ 2015. Shear-mediated contributions to the effective properties of soft acoustic metamaterials including negative index. Sci. Rep. 5, 18562 (doi:10.1038/srep18562) [PMC free article] [PubMed]
19. Norris AN. 1985. A differential scheme for the effective moduli of composites. Mech. Mat. 4, 1–16. (doi:10.1016/0167-6636(85)90002-X)
20. Haberman MR, Berthelot YH, Jarzynki J, Cherkaoui M 2002. Micromechanical modeling of viscoelastic voided composites in the low-frequency approximation. J. Acoust. Soc. Am. 112, 1937–1943. (doi:10.1121/1.1509424) [PubMed]
21. Christensen RM. 1982. Theory of viscoelasticity, 2nd edn New York, NY: Dover Publications, Inc.
22. Dvorak G. 2013. Micromechanics of composite materials. Berlin, Germany: Springer.
23. Haberman MR, Berthelot YH, Cherkaoui M 2006. Micromechanical modeling of particulate composites for damping of acoustic waves. J. Eng. Mat. Tech. 128, 320–329. (doi:10.1115/1.2204943)
24. Hill R. 1965. A self-consistent mechanics of composite materials. J. Mech. Phys. Solids 13, 213–222. (doi:10.1016/0022-5096(65)90010-4)
25. Hershey AV. 1954. The elasticity of an isotropic aggregate of anisotropic cubic crystals. ASME J. Appl. Mech. 21, 236–240.
26. Kröner E. 1958. Berechnung der elastishen konstanten des vielkristalles aus den konstanten des einkristalls. Z. Phys. 151, 504–518. (doi:10.1007/BF01337948)
27. Fang N, Xi D, Xu J, Ambati M, Srituravanich W, Sun C, Zhang X 2006. Ultrasonic metamaterials with negative modulus. Nat. Mat. 5, 452–456. (doi:10.1038/nmat1644) [PubMed]
28. Fassi-Fehri O. 1985. Le problème de la paire d'inclusions plastiques et hétérogènes dans une matrice anisotrope. Metz, France: Université de Metz.
29. Haberman MR. 2007. Design of high loss viscoelastic composites through micromechanical modeling and decision based materials design. Atlanta, GA: Georgia Institute of Technology.

Articles from Proceedings. Mathematical, Physical, and Engineering Sciences are provided here courtesy of The Royal Society