|Home | About | Journals | Submit | Contact Us | Français|
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.
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  for an introductory discussion) by using dynamic Eshelby tensors . 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 . In another micromechanical approach, Gonnella & Ruzzene  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 . For example, Baird et al.  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.  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  for a review of some models, and Lakes  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 . While Sabina & Willis  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’ , 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 . 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.
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.
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 , 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 .
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.
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.
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. ). Similarly, two dots sum over the last two indexes of the left tensor and the first two indexes of the right tensor (e.g. ). 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 is the momentum density. The operator diviis the divergence over the ith dimension (e.g. ), 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, is the linear strain tensor, ρ is the density, is the displacement field, and is the displacement gradient field (i.e. ). The raised ‘T’ denotes a transpose which switches the order of the first two dimensions of the tensor (e.g. ). 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  should be employed. The constitutive relations for the stress and momentum fields in a Willis material are
where and 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 . Note that the Willis equations reduce to the standard constitutive equations for 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 , 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 is the frequency domain Green's function and is the three-dimensional Dirac delta function. The Green's function gives the displacement field at the position due to a point force source at the origin (or at a source location if the arguments of and are replaced by ). Furthermore, is material, geometry and frequency dependent, in general.
The displacement field may then be implicitly written as
where 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
Thus, the integral equation for the velocity may be written as
Another property of convolutions, , 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. ).
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 which is equal to zero if is in the matrix and equal to one if is in the inhomogeneity. Now define ΔX such that for , 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′2is the divergence of the second dimension with respect to . Note that the solution to the homogeneous boundary-value problem is still accounted for when . 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 .
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, kMa1, 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 . 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 and , where is the characteristic specific acoustic impedance of the material and || denotes the characteristic amplitude of the argument. The interpretation of and becomes apparent when considering the order of magnitude of the terms in the constitutive equations (2.3):
where k is the characteristic wavenumber. Thus, and represent the relative importance of Willis coupling in the constitutive equations. By reciprocity and passivity, we also know that [15,16]. If we define and H=G|C|a, we may write equation (2.8) in the dimensionless form
where partial differentiation is now with respect to and we have assumed , which is equivalent to assuming . Note that in the case of strong coupling, or , the wavenumber is approximately . The characteristic amplitudes of the three coefficients in equation (2.17) are, in order from left to right, 1, , 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 . 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 X0 is the homogeneous field X0 averaged over the domain of the inclusion and
The expressions for and 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 , then we may write and . Also, note that |T1|~1/|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 kMaWM1 and (kMa)21 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 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 IIijkl=(δikδjl+δilδjk)/2 is the fourth-order symmetric identity tensor.
Because the tensors Aε, , 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, calculates the microscopic strain due to the macroscopic velocity, 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.
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 , where the superscript eff denotes the effective, or macroscopic, property and ϕΩI/[ΩI+ΩM] is the inhomogeneity volume fraction. The averaged fields YM and YI 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, ε and in equation (2.25) may be used to represent Y I. This is the dilute approximation, which is expressed mathematically as kMakML, 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/31.
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, . This is equivalent to the classical Eshelby method and is valid up to O(ϕ) . Thus,
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 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 and Δρ=0, which is a slight modification from the original observation made by Willis , 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].
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 , the DEM equations which extend the present quasi-static dilute effective material approximation may be written as
with the conditions Ceff(ϕ=0)=CM, , 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 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 . This inequality, which violates reciprocity [15,16], is a consequence of the quasi-static approximation. Despite this fact, the prediction of 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 , 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 . 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.
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 , 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.
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=1GPa, bulk modulus KM=3GPa, and density ρM=700kgm−3 containing 0.1mm radius spherical inhomogeneities with shear modulus μI=10GPa, bulk modulus KI=30GPa and density ρI=2000kgm−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 50Hz.
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 , 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.
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 50Hz, the inhomogeneity radius is 0.1mm, 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.
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.
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 50Hz.
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.
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 , 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 . The choice of F(ω) follows the form given by Fang et al.  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 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 , 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 T:ΔC≈I. 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.
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 (kMa1, 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(kMa), only the strain localization is accurately predicted, and so only one of the Willis coupling tensors, , 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.
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 , define the three-dimensional spatial Fourier transform (at a single frequency) with the pair
where ki=kχi is the wavevector, k is the wavenumber, and χi is the unit direction-cosine vector. Note the transforms have the property . For later use, note
(a) Derivation of for spherical inhomogeneities
The Fourier transform of the quasi-static Green's function equation given in equation (2.18) yields , which may be solved for Green's function as , where Pik=χlχjCijkl.
The tensor T1 as defined in equation (2.21a) may be expanded as , where
Ω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
which must be integrated numerically, the tensor may be written as
(b) Derivation of for spherical inhomogeneities
The tensor as defined in equation (2.21b) may be written in index notation as
Substituting Gik,j for its Fourier transform and performing similar integrals as in the analysis of , we may write , where
Thus, . Note that since , the value of must be zero for centre-symmetric geometries, and so .
(c) Extension to ellipsoidal inhomogeneities
Consider an ellipsoidal inhomogeneity defined such that the outer boundary of the inhomogeneity is expressed as
Define a new coordinate system Xi such that xi=ΦijXj, where
This suggests Ki should be defined such that . In this coordinate system, the inhomogeneity is a sphere of radius a1 and xiki=XiKi. The Jacobian determinant of the spatial transformation is , and the Jacobian determinant of the wavevector transformation is . Then may be written as
where Ωk is the k-space volume. Recall χi is the vector of direction cosines of ki. Thus, , where χi′ is the vector of direction cosines of Ki. Defining P′ikχ′lχ′jC′ijkl, where , such that , results in
which becomes , where
Following similar steps, , 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 and . 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
and the inverse may be written as . Note is a symmetric tensor. The needed integrals are then written
Then for an isotropic matrix with spherical inhomogeneities and
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.
We have no competing interests.
This work was supported by ONR through MURI grant no. N00014-13-1-0631.