## Abstract

This article presents a method for the homogenization of three-dimensional periodic elastic composites. It allows for the evaluation of the averaged overall frequency-dependent dynamic material constitutive tensors relating the averaged dynamic field variable tensors of velocity, strain, stress and linear momentum. Although the form of the dynamic constitutive relation for three-dimensional elastodynamic wave propagation has been known, this is the first time that explicit calculations of the effective parameters (for three dimensions) are presented. We show that for three-dimensional periodic composites, the overall compliance (stiffness) tensor, as produced directly by our formulation, is Hermitian, regardless of whether the corresponding unit cell is geometrically or materially symmetric. Overall, mass density is shown to be a tensor and, like the overall compliance tensor, always Hermitian. The average strain and linear momentum tensors are, however, coupled, and the coupling tensors are shown to be each others' Hermitian transpose. Finally, we present a numerical example of a three-dimensional periodic composite composed of elastic cubes periodically distributed in an elastic matrix. The presented results corroborate the predictions of the theoretical treatment illustrating the frequency dependence of the constitutive parameters. We also show that the effective properties calculated in this paper satisfy the dispersion relation of the composite.

## 1. Introduction

Recent interest in the character of the overall dynamic properties of composites with tailored microstructure necessitates a systematic homogenization procedure to express the dynamic response of an elastic composite in terms of its average effective compliance and density. The elastostatic response of composites has been long understood to be non-local in space (Beran 1968; Willis 1983; Diener *et al.* 1984; Bakhvalov & Panasenko 1989), but in the context of inhomogeneous elastodynamics, the effective constitutive relations are non-local in both space and time (Willis 1981*a*,*b*). Homogenization for calculating these overall dynamic properties of composites, based on the integration of the field variables, has been proposed by a number of researchers. For electromagnetic waves, see, for example, Smith & Pendry (2006), Amirkhizi & Nemat-Nasser (2008), Bensoussan *et al.* (1978) and Sihvola (1999). For elastodynamic waves, Willis (2009) has presented a homogenization method based on an ensemble averaging technique of the ‘Bloch’ reduced form of the wave propagating in a periodic composite; see also Shuvalov *et al.* (2011) and Willis (2011).

In the present paper, we propose a method of homogenization of three-dimensional microstructured elastic composites, which directly provides the overall frequency-dependent dynamic material parameters. This method is the three-dimensional generalization of the one-dimensional homogenization scheme presented by Nemat-Nasser & Srivastava (2011). The one-dimensional homogenization scheme was shown to be equivalent to the field variable integration scheme of Nemat-Nasser *et al.* (2011) in the limit of Bloch wave propagation.

Our formulation is such that neither the pointwise values of the field variables nor a Green function are required for the calculation of the effective properties. We replace the heterogeneous unit cell by a homogeneous one with a conveniently chosen uniform elasticity tensor and mass density (the final results are independent of this choice). Then, we introduce distributed eigenstrains and eigenmomentums, or eigenstresses and eigenvelocities, such that the pointwise values of the field variables remain the same as they would be in the original heterogeneous solid. However, our formulation uses only the *volume averages* of these eigenfields to calculate the overall constitutive tensors. Unlike pointwise values, these volume averages can be accurately estimated using crude approximations.

If desired, the dispersion relation of the composite can be extracted from the final homogenized constitutive parameters. For this, we constrain the overall response to follow a Bloch-form wave motion, leading to an eigenvalue problem in terms of the overall frequency-dependent constitutive parameters.

In the following, we outline our homogenization approach for a general three-dimensional periodic elastic composite. The resulting averaged dynamic constitutive parameters are tensorial in nature and the average strain tensor is coupled with the average momentum tensor. Such a form of the averaged constitutive relation, where the constitutive parameters (including mass density) are tensors and where the average strain (average stress) is coupled with average linear momentum, has been predicted in the literature (Willis 1997; Milton & Willis 2007; Willis 2009; Shuvalov *et al.* 2011; and references therein), but nowhere in the literature have these effective properties been explicitly calculated for full three-dimensional wave propagation or vibration. In general, the effective properties for the dynamic problem are not uniquely determined, but they become unique in the presence of incompatible strains (Willis 2011). Our method (which uses a comparison medium) directly yields these ‘unique’ parameters (J. R. Willis 2010, personal communication). This paper is a detailed exposition on calculating these ‘unique’ effective parameters from the comparison medium. Furthermore, we prove that the coupling parameters arising in the averaged constitutive relations given by our method are the Hermitian transpose of each other, regardless of the material or geometric asymmetries.

Finally, it should be pointed out that our approach can also be used to calculate the pointwise values of the field variables, although this has no bearing on the calculation of the effective properties by our method.

## 2. Microstructural homogenization of periodic composites

Here, we present a homogenization method based on micromechanical consideration of the volume averages of the field variables, viewed as measurable macroscopic physical quantities. The treatment presented here is the dynamic equivalent of the static homogenization theory already firmly established in the literature (see Nemat-Nasser & Hori 1999; Torquato & Haslach 2002 and references therein for a comprehensive account). We express the solution to the elastodynamic equations of motion as the sum of the volume average and an additional deviation field owing to the heterogeneous composition of the unit cell, 2.1 The aim is to derive a set of constitutive relations for the overall averaged parts of the field variables, using the local elastodynamic equations of motion and constitutive relations. This then provides the homogenized frequency-dependent material parameters.

Consider harmonic waves in an unbounded elastic composite consisting of a collection of bonded, identical unit cells (*Ω*={*x*_{i}:−*a*_{i}/2≤*x*_{i}<*a*_{i}/2; *i*=1,2,3}), which repeat themselves in all directions, and hence constitute a periodic structure. In view of the periodicity of the composite, we have *ρ*(**x**)=*ρ*(**x**+*m*′**I**_{β}) and **C**(**x**)=**C**(**x**+*m*′**I**_{β}); here, *m*′ is an integer, *ρ*(**x**) is the density, **C**(**x**) is the fourth-order tensor of the modulus of elasticity whose inverse is the compliance tensor **D**(**x**) and **I**_{β},*β*=1,2,3, denote the three vectors that form a parallelepiped enclosing the periodic unit cell. For time harmonic waves with frequency *ω*, the field quantities are proportional to e^{±iωt}. For waves with wavevector **q**=*q*_{i}**e**_{i}, where **e**_{i} is the unit vector in the *i*th direction and the Einstein summation convention applies, the Bloch representation of the field variables takes the following form:
2.2
where represents the field variables, stress (), strain (), momentum () or velocity (), whereas **Q** represents their periodic parts (** σ**,

**,**

*ε***p**, ). The representation, equation (2.2), separates the time harmonic and macroscopic factor from the microscopic part of the field variables. We emphasize here that the frequency and the wavevector,

*ω*and

**q**, are, at this point, unrelated and arbitrary. Indeed, with

**q**equal to zero and

*ω*prescribed, the corresponding dynamic effective properties are obtained.

The local conservation and kinematic relations are
2.3
where . The corresponding local constitutive relations are
2.4
where **D**(**x**)=**C**(**x**)^{−1} is the tensor of compliance and *ρ*(**x**) is the density of the material. These local material parameters represent the structure and composition of the unit cell.

Now, we replace the heterogeneous unit cell with a homogeneous one having a suitable positive uniform density *ρ*^{0} and a positive-definite compliance **D**^{0}=[**C**^{0}]^{−1} with the usual symmetries. These reference material parameters can be chosen for convenience without affecting the final overall average properties. In order to reproduce the strain and momentum of the actual unit cell, field variables eigenstrain, **E**(**x**), and eigenmomentum, **P**(**x**), are introduced. These quantities are then calculated, using the basic local field equations and constitutive relations. The idea stems from the polarization stress or strain that was originally proposed by Hashin (1959), and further developed by Hashin & Shtrikman (1962*a*,*b*), Hashin (1963) and later by others, in order to construct energy-based bounds for the composite's overall elastic moduli. The basic tool in these works has been the result, obtained by Eshelby (1957) in three dimensions and earlier by Hardiman (1954) in two dimensions, that the stress and strain are constant within an ellipsoidal (elliptical in two dimensions) region of an infinitely extended uniform elastic medium when that region undergoes a uniform transformation corresponding to a uniform inelastic strain.

Here, however, we present a different tool that can also be used to calculate the pointwise values of the elastodynamic field variables to any desired degree of accuracy; our averaging method, however, requires only volume averages. For this, we require that the actual values of the field variables at every point within the homogenized and the original heterogeneous unit cell be exactly the same. To ensure this, we require that the following *consistency conditions* hold at every point within the unit cell (**x** dependence implicit in all field variables and eigenfields):
2.5
The eigenstrain and eigenmomentum fields are zero in regions where the material properties of the heterogeneous unit cell are equal to the chosen uniform material properties **D**^{0} and *ρ*^{0}. From equations (2.3) and (2.5), we have
2.6
and
2.7
Since the stress and displacement fields (**Q**) are periodic (also for representative volume element (RVE), viewed as a unit cell), they can be expanded in a spatial Fourier series,
2.8
2.9
2.10
2.11
and
2.12
where Greek indices are not summed. In the above equations, 〈**Q**〉 represents the averaged value of the field variable over the unit cell and appears in its macroscopic description, and **Q**^{p}**(x)**, which is periodic with zero mean value, represents the local deviations from the average value. Equations (2.6) and (2.7) become
2.13
and
2.14
where ** ζ**=

**+**

*ξ***q**. For the case of an isotropic reference material, we have 2.15 where

*λ*

^{0}and

*μ*

^{0}are the Lamé constants. Using the isotropic stiffness tensor, the tensors in equations (2.13) and (2.14) can be inverted to give 2.16 where explicit expressions for tensors

**,**

*Φ***,**

*Θ***and**

*Ψ***are given in appendix A. Now, the stress and velocity fields can be expressed as a sum of their average and zero mean periodic components, 2.17 and 2.18 where and 〈**

*Γ***〉 are the average values of the velocity and stress fields, respectively, taken over a unit cell.**

*σ*To make the homogenized unit cell pointwise equivalent to the original heterogeneous unit cell, the homogenizing fields are required to satisfy the following consistency conditions:
2.19
and
2.20
where **D** and **D**^{0} are the compliance tensors of the actual and the reference materials, respectively. The periodic parts of the velocity and stress fields, from equations (2.17) and (2.18), are now substituted into the above equations. This gives a set of two coupled integral equations which yields the required homogenizing stress and velocity fields that exactly and fully replace the heterogeneity in the original medium.

To obtain the overall effective properties, we only need to calculate the average quantities. To this end, average equations (2.19) and (2.20) over a unit cell (or an RVE) and obtain
2.21
and
2.22
To calculate 〈**E**〉 and 〈**P**〉, we divide the unit cell into subregions, *Ω*_{α}. Then, we average the periodic fields over each such subregion to obtain
2.23
2.24
2.25
We now replace the integrals in equations (2.23) and (2.24) by their equivalent finite sums and set
2.26
Equations (2.23) and (2.24) then yield the following expressions:
2.27
and
2.28
where the repeated index, *β*, is summed over the number of subregions, , and Greek indices serve as labels for tensors rather than components of a particular tensor. The coefficient tensors in the above equations are defined by
2.29
In these equations, *α* and *β* are not summed. Averaging the consistency conditions over each subregion *α* and using equations (2.27) and (2.28), we have
2.30
Equations (2.27) and (2.28) can now be inverted to express the eigenstrain tensors **E**^{β} and eigenmomentum tensors **P**^{β} in terms of the average stress 〈** σ**〉 and average velocity tensors. In addition, these piecewise constant eigenfields can be used to express the unit-cell-averaged eigenfields in terms of the average stress and average velocity tensors. The solution is formally written as
2.31
It should be stressed here that only the averages of the eigenfields are required for the determination of the effective properties. This, in turn, means that unlike the field integration method of Nemat-Nasser

*et al.*(2011) or Smith & Pendry (2006), the current formulation does not employ the pointwise elastodynamic solution. The formulation, however, is consistent with the elastodynamic problem, and the dispersion relation and the eigenvectors of the composite can be extracted from it (§2.1). The averaged consistency conditions (2.21) and (2.22) are now expressed as 2.32 and 2.33

Equations (2.32) and (2.33) are our final constitutive relations for the homogenized composite. It is shown in appendix B that the effective compliance tensor is such that and the effective density tensor is such that . It is also shown that the coupling tensors have the relation .

The energy corresponding to the averaged field quantities is given by 2.34 and 2.35 and is real-valued. Taking into account the symmetries of the homogenized constitutive parameters, note that the averaged constitutive relations of equations (2.32) and (2.33) can have 45 independent constants at the maximum. Depending on the material properties of the constituents and the symmetry of the unit cell, the number of independent constants will vary. The essential relations of the constitutive tensors proved in appendix B will hold, regardless of the material anisotropy or directional asymmetry.

Instead of considering the micromechanical formulation in terms of eigenstrains and eigenmomentums, we may instead consider it in terms of eigenstresses and eigenvelocities. This would lead to an effective constitutive relation that can be expressed as (Willis 2009, 2011) 2.36 and 2.37 For Bloch-wave propagation, the two approaches are equivalent and the above constitutive parameters can be directly extracted from equations (2.32) and (2.33). Specifically, , , and . For general boundary conditions, however, the two homogenization schemes may not be equivalent. Willis has shown that the constitutive relations expressed in the above form have a self-adjoint structure, i.e. . Using the properties of the effective parameters proved in appendix B, it can be shown that the self-adjointness of the structure of the constitutive relations of equations (2.36) and (2.37) emerges identically from our formulation.

### (a) Bloch-form overall spatial variation and dispersion relations

The wavenumber **q** and frequency *ω* have been assumed to be independent until now. We now consider a special case of an infinite homogenized elastic solid with a periodic microstructure, and seek conditions under which it supports periodic waves of the following Bloch-form spatial variation:
2.38
The averaged field equations are given by the kinematical relations (2.3) as
2.39
These equations are combined with the constitutive relations (2.21) and (2.22) to eliminate three out of the four field variables and arrive at a single equation of the following form:
2.40
where **R** is a second-order tensor whose components are functions of the wavenumber **q** and the frequency *ω*, as well as the geometrical and material properties of the unit cell, which are reflected through the overall constitutive relations (2.32) and (2.33). For non-trivial solutions characterized by a velocity vector with non-zero components, the above equation is satisfied for specific **q**–*ω* pairs. These **q**–*ω* pairs constitute the dispersion relation of the composite for **q**-propagation. The associated non-trivial velocity vector (eigenvector) indicates the polarization of the Bloch wave. In this way, the formulation presented in this paper can be used to independently calculate the dispersion relation of the composite.

In §3, we calculate the dispersion relation of a specific three-dimensional periodic composite and compare the results with published results in the literature. We then present effective property calculations for two points on the dispersion curve, thereby explicitly demonstrating both the structure and the frequency dependence of the constitutive relation.

## 3. Numerical example: elastic cuboids periodically distributed in an elastic matrix

We consider a composite with a three-dimensional periodic microstructure composed of the periodically distributed unit cell shown in figure 1. The wave is propagating in the *x*_{1} direction so that **q**=*q***e**_{1}. The Poisson ratio of the matrix *ν*_{M}=0.35 and the inclusion *ν*_{I}=0.3. The ratio of the elastic modulus *E*_{I}/*E*_{M}=50 (*E*_{M}=1 GPa) and the ratio of the density *ρ*_{I}/*ρ*_{M}=3 (*ρ*_{M}=1000 kg m^{−3}). The inclusion is a cuboid of sides *b*_{i} in the directions *x*_{i} and the matrix is a cuboid of sides *a*_{i} in the directions *x*_{i}. For the specific problem of figure 1, we have assumed *a*_{1}/*a*_{2}=*a*_{1}/*a*_{3}=1, *b*_{1}/*a*_{1}=0.8, *b*_{2}/*a*_{2}=0.9 and *b*_{3}/*a*_{3}=0.7.

We now compare dispersion calculations from two different methods for the first shear vertical (displacement primarily in the *x*_{2} direction) and shear horizontal (displacement primarily in the *x*_{3} direction) modes. Table 1 gives the normalized frequency calculations (for different values of normalized wavenumber) from both the mixed-variational formulation and the micromechanical method given in this paper (see Minagawa & Nemat-Nasser (1976) for definitions of normalized wavenumber and normalized frequency).

The mixed-variational formulation is a convergent and accurate method for determining the dispersion relation of periodic composites (see Babuska & Osborn (1978). A brief summary is provided in appendix C). The results presented in table 1 under the mixed-variational calculations match previously published results in Minagawa & Nemat-Nasser (1976). The corresponding micromechanical calculations are carried out by discretizing the cuboidal inclusion into eight smaller cuboids. A total of 15^{3}−1 Fourier terms were used in the expansion of the periodic parts of the field variables. The results show good correspondence with the independently calculated results from the mixed-variational formulation, and indicate that the effective parameters calculated from the micromechanical method satisfy the dispersion relation of the composite.

Tables 2 and 3 show the effective properties calculated from the micromechanical homogenization method at two different frequency–wavenumber pairs for the shear vertical mode. As has been shown above, these dynamic effective properties satisfy the dispersion relation of the composite. It should be noted that all the components of the four effective tensors are real. This is due to the fact that the unit cell is symmetric in all three orthogonal directions. The essential relation between the coupling tensors () is verified from the present calculation. The effective compliance tensor has both major and minor symmetries () owing to the symmetry of the unit cell with respect to the coordinate directions. If the unit cell is asymmetric, then, in general, we should have . The effective density matrix is diagonal and indicates that there is no coupling between the average momentum in a particular direction and the average velocity in an orthogonal direction. Tables 2 and 3 also serve to clarify the frequency dependence of the effective parameters.

## 4. Conclusions

A method for the homogenization of three-dimensional periodic elastic composites is presented. It allows for the evaluation of the overall dynamic material constitutive tensors relating the averaged dynamic field variable tensors of velocity, strain, stress and linear momentum. The method does not require the *a priori* calculation of the solution to the elastodynamic problem nor the pointwise values of any field variable. To calculate the effective constitutive tensors, only the volume averages of the homogenizing eigenfields are required. These volume averages can then be accurately estimated using crude approximations. When the solution is constrained to follow a Bloch-form wave motion, then the dispersion relations of the composite can be extracted from the final effective properties by solving an eigenvalue problem. The coupled form of the constitutive relation as proposed by Milton & Willis (2007) emerges naturally from the present method. We have also shown that the matrices corresponding to the effective compliance and density tensors that are produced by our method are Hermitian and that the coupling tensors are Hermitian transposes of each other. Finally, we have presented a numerical example showing that the effective properties calculated from the present method satisfy the dispersion relation of the composite. The example also serves to clarify the structure, symmetries and frequency dependence of the constitutive parameters. In addition (and unrelated to the averaged properties), the current approach can be used to calculate the pointwise values of the field variables.

## Acknowledgements

The authors are grateful to Prof. John R. Willis for valuable comments. This research has been conducted at the Center of Excellence for Advanced Materials (CEAM) at the University of California, San Diego, under DARPA AFOSR Grants FA9550-09-1-0709 and RDECOM W91CRB-10-1-0006 to the University of California, San Diego.

## Appendix A. Explicit expressions for tensors

Equation (2.16)_{1} in index notation is
A1
where ** Φ** and

**for the isotropic reference matrix are given by A2 and A3 In the above expressions, is the longitudinal wave velocity and is the shear wave velocity.**

*Θ*Equation (2.16)_{2} in index notation is
A4
where
A5
and
A6
where
It can be seen from the symmetries of the **S** tensor that the above tensors satisfy the following properties:
A7

## Appendix B. Mathematical structure of the constitutive relation

For each tensor **M**^{αβ} in equations (2.27) and (2.28), we have
B1
This property also holds for tensors and . To facilitate inversion and solution, equations (2.30) are now expressed in their equivalent matrix form,
B2
In the above equations, [** Γ**] is a square matrix and [

**] is a square matrix. These matrices are composed of the components of the tensors and appropriately placed at the matrix locations.**

*Φ*It can be shown that for any *α* and *β*, we have
B3
where , for instance, is the *pqrs* component of the tensor . The above result signifies that the matrices [** Γ**] and [

**] are Hermitian. We also note the following identity for the matrices [**

*Φ***] and [**

*Ψ***]: B4 More generally, denoting a Hermitian transpose by †, we have the following identities for the matrices: B5**

*Θ**To expedite notation*, we omit the square brackets denoting a matrix in further calculations.

*It must be stressed*that while unbracketed quantities in the main text represent tensors, they represent matrices in this appendix. This is done purely for the representational ease of the essential proofs that follow. We express the solution to equation (B2) in the following form (matrix form of equation 2.31): B6 where B7

Equation (B6) expresses the vectors of eigenstrain and eigenmomentum in each subdivision in terms of the vector of average stress and average velocity. It can be seen that the vectors of the average quantities in equation (B6) are composed of times repeated vectors of the averaged quantities. Therefore, equation (B6) can be condensed to express the vectors of eigenstrain and eigenmomentum in each subdivision in terms of a 6×1 vector of the average stress and a 3×1 vector of the average velocity. Furthermore, we can also average the vectors of eigenstrain and eigenmomentum over all the subregions to finally get B8 where B9 where and

Now, we average the consistency conditions over the entire unit cell to express the average strain and average momentum in terms of the average stress and average velocity tensors. Noting that the average of the periodic parts vanishes when taken over the entire unit cell, we have B10 and B11 The effective parameters are given by B12 The matrix representation of the effective constitutive relation (B10) and (B11) can be transformed into the tensorial representation of the main text (2.32) and (2.33) without any loss of generality.

(*a*) *Properties of the constitutive parameters*

Considering equation (B5), it can be shown that 〈** Λ**〉 and 〈

**〉 are Hermitian. This implies that and are also Hermitian.**

*Ω*Taking a Hermitian transpose of equation (B12^{3}),
B13
proving that . Since the matrices and are Hermitian, we have and , where ‘*’ denotes a complex conjugate. This means that the scalar given by has a complex conjugate . Since the pairs *i*,*j* and *k*,*l* are symbols upon which summation is carried out, we find that *a*=*a** or that *a* is real. Similarly, can be shown to be a real scalar. The relation implies that .

## Appendix C. Brief overview of the mixed-variational formulation

Consider harmonic waves in an unbounded periodic elastic composite consisting of a collection of unit cells, *Ω*. In view of periodicity, we have *ρ*(**x**)=*ρ*(**x**+*m*′*I*^{β}) and *C*_{jkmn}(**x**)=*C*_{jkmn}(**x**+*m*′*I*^{β}), where *m*′ is any integer and *I*^{β},*β*=1,2,3, denote the three vectors that form a parallelepiped enclosing the periodic unit cell.

For time harmonic waves with frequency *ω* (*λ*=*ω*^{2}), the field quantities are proportional to e^{±iωt}. The field equations become
C1

For harmonic waves with wavevector **q**, the Bloch boundary conditions take the form
C2
for **x** on ∂*Ω*, where **t** is the traction vector.

To find an approximate solution of the field equations (C1) subject to the boundary conditions (C2), we consider the following expressions:
C3
and
C4
where the approximating functions *f*^{(αβγ)} are continuous and continuously differentiable, satisfying the Bloch periodicity conditions. As shown in Minagawa & Nemat-Nasser (1976), the eigenvalues are obtained by rendering the following functional stationary:
C5
where , with ‘*’ denoting complex conjugate, and *D*_{jkmn} are the components of the elastic compliance tensor, the inverse of the elasticity tensor *C*_{jkmn}.

Substituting equations (C3) and (C4) into equation (C5) and equating to zero the derivatives of *λ*_{N} with respect to the unknown coefficients and , we arrive at the following set of linear homogeneous equations:
C6
There are (*M*_{p}=2*M*+1) equations in equation (C6)_{2} for a general three-directionally periodic composite. They may be solved for in terms of and the result substituted into equation (C6)_{1}. This leads to a system of linear equations. The roots of the determinant of these equations give estimates of the first eigenvalue frequencies. The corresponding eigenvectors are , from which the displacement field within the unit cell is reconstituted.

- Received July 22, 2011.
- Accepted August 26, 2011.

- This journal is © 2011 The Royal Society