## Abstract

A new theory of elastodynamic homogenization is proposed, which exploits the integral equation form of Navier's equations and relationships between length scales within composite media. The scheme is introduced by focusing on its leading-order approximation for orthotropic, periodic fibre-reinforced media where fibres have arbitrary cross-sectional shape. The methodology is general but here it is shown for horizontally polarized shear (SH) wave propagation for ease of exposition. The resulting effective properties are shown to possess rich structure in that four terms account separately for the physical detail of the composite (associated with fibre cross-sectional shape, elastic properties, lattice geometry and volume fraction). In particular, the appropriate component of Eshelby's tensor arises naturally in order to deal with the shape of the fibre cross section. Results are plotted for circular fibres and compared with extant methods, including the method of asymptotic homogenization. The leading-order scheme is shown to be in excellent agreement even for relatively high volume fractions.

## 1. Introduction

The theory of elastodynamic homogenization has a rich history, and many methods have been proposed in order to find the effective properties of composite media in the true homogenization regime when the wavelength of the propagating wave *λ* is much larger than the characteristic length scale of the microstructure *q*. Thus, on defining , the leading-order behaviour in *ϵ* is static and higher order terms (in *ϵ*) correspond to frequency dependence and dispersive effects.

For periodic media, the statics problem has been well studied by the method of asymptotic homogenization (MAH) (see Bakhvalov & Panasenko 1989; Parton & Kudryavtsev 1993). In a series of papers, Sabina *et al*. (2002 and references therein) employed Weierstrassian elliptic functions to solve static problems using the MAH for materials of various symmetry and for transversely isotropic phases. Eigenstrain methods (see Nemat-Nasser & Hori 1999) have also been applied in the static context by Jiang *et al.* (2004).

Extensions to the multipole scheme developed in the pioneering work of Rayleigh (1892) on the effective electrical properties of inhomogeneous media have been applied in many contexts since then. In particular, it was employed first in linear elasticity by McPhedran & Movchan (1994) and in an elastodynamic context by Zalipaev *et al.* (2002) in order to find stop bands of periodic fibre-reinforced elastic structures. In the low-frequency limit, the scheme also provides predictions of the effective moduli, although it can become highly unstable in the large volume fraction regime. The MAH may also be applied in a low-frequency elastodynamic context as was shown by Parnell & Abrahams (2006), who determined the effective quasi-static longitudinal shear properties of periodic monoclinic fibre-reinforced media. A particularly attractive feature of this scheme is the fact that an anisotropic wave equation emerges naturally, thereby implying the elastic symmetry of the composite in question. Recently, Andrianov *et al.* (2008) have analysed the dispersive corrections to the classical asymptotic results.

For random media, several self-consistent schemes have been proposed (e.g. see the work by Sabina & Willis 1988). These schemes are often difficult to justify mathematically or physically although steps to address this have recently been taken by Markov (2001) and Kanaun & Levin (2003). The method of multiple scattering initiated by Foldy (1945) and developed by Waterman & Truell (1961) and Fikioris & Waterman (1964) has also been successfully applied in elastodynamic and acoustic contexts (see Bose & Mal 1974 and Linton & Martin 2005, respectively) and by employing the T-matrix scattering theory by Varadan *et al.* (1978). The multiple-scattering formulation allows correlation functions to be introduced into the analysis, although it requires a closure condition (the quasi-crystalline approximation being the usual choice), which itself is difficult to justify rigorously. The quasi-static limit of the multiple-scattering prediction in the *well-stirred* limit (see Varadan *et al.* 1978) corresponds exactly to the special *composite cylinders* result of Hashin & Rosen (1964), when a particular limit corresponding to *no gaps between composite cylinders*, is taken. However, although multiple-scattering theory is well suited to isotropic phases, it becomes extremely complicated for anisotropic phases. In that respect, as noted by Willis (1980), an integral equation approach is preferable. In a static context, successful calculations of effective moduli for composite media with complex microstructure via numerical solutions of integral equations have been calculated by Greengard & Helsing (1998), albeit with various assumptions on the resulting effective moduli.

Recently, a new integral equation method (IEM) of homogenization was proposed by Brazier-Smith (2002), but this method was shown to be valid only for *dilute dispersions* by Parnell (2004). This scheme will henceforth be called the dilute dispersion integral equation (DDIE) method. This approach has, however, motivated a new integral equation formulation for homogenization in the non-dilute volume fraction regime, which we introduce here in the context of fibre-reinforced media where fibres are arranged on a *periodic* lattice. This method presents a number of advantages over the existing homogenization schemes discussed above.

For the purposes of introducing this method, we consider the simplest elastic wave propagation problem: the scalar problem of horizontally polarized shear (SH) wave propagation through the fibre-reinforced composite, which allows the determination of longitudinal shear properties of the medium. Furthermore, we restrict attention to those lattice geometries that result in (at most) orthotropic elastic symmetry on the macroscale. We assume that each fibre within the composite is identical and that all phases that make up the composite are isotropic. We stress, however, that the case of families of fibres with different cross sections and arbitrary (at most monoclinic) anisotropy can also be dealt with by this method.

In this paper, we introduce the integral equation formulation by considering its leading-order system (in volume fraction) alone; note, however, that this is *not* equivalent to a dilute volume fraction approximation since interaction effects are incorporated. In particular, in contrast to dilute approximations, we do *not* assume that the wavelength of the propagating wave is less than the characteristic length scale of the separation of fibres. We show that the leading-order result is itself an extremely good approximation to the effective properties, showing an excellent agreement with extant methods even for fairly high volume fractions. We discuss how the approach can be applied to materials having fibres of arbitrary cross section and plot results for the common case of circular cross sections.

In §2 we introduce appropriate notation and the governing equations and in §3 the integral equation methodology is presented. Displacement and displacement-gradient *moments* (weighted averages of field variables) are defined and the governing equations for the leading-order moments are found (this is the so-called leading-order system referred to above). Wave-like ansatz are posed for the moments, thereby introducing the effective wavenumber *γ*. A separate term that accounts completely for the shape of the fibre cross section arises and is directly related to Eshelby's tensor. The leading-order (homogeneous) system is given in §3*c* and the condition allowing us to determine the effective wavenumber follows. We note, in particular, that this expression incorporates terms associated with volume fraction, elastic properties, lattice summations and a shape tensor.

In §4, we show how the lattice summations may be determined when fibres are positioned on rectangular (orthotropic) and hexagonal (transversely isotropic) lattices and determine the form of the shape tensor explicitly. The effective quasi-static properties of the composite are determined in §5; results are given in the most general case and we plot the solutions in the case when the fibre cross section is circular, comparing these with extant methods. We show that, even at leading order, the proposed integral equation scheme is reliable and provides a new algorithmic approach to elastodynamic homogenization. We close in §6 with conclusions.

## 2. Governing integral equations

We shall work in Cartesian coordinates where the *x*_{3} axis runs parallel to the fibres and time-harmonic SH waves propagate in the *x*_{1}*x*_{2} plane, polarized along the *x*_{3} axis with circular frequency *ω*. We define the lattice vector(2.1)which is the position vector of the centre of the th periodic cell of the material (). For ease of exposition here, we shall consider the case where the periodic cell contains only one fibre and thus ** R** also denotes the

*origin*of the th fibre (i.e. some convenient point chosen inside the fibre domain). We see from (2.1) that

*q*is the characteristic length scale of the periodic cell (and thus the microstructure) and are vectors whose components are . Again for ease of exposition, we shall consider distributions of fibres which would result in an orthotropic material. All geometries of this type can be described by lattice vectors of the form1(2.2)with , and therefore the cell is a periodic parallelogram (figure 1) of area .

As shown in figure 1, the th periodic cell comprises isotropic host phase of mass density and shear modulus (and thus wavenumber where is the square of the host phase speed) and a fibre of general cross section (not restricted to being circular) having domain , mass density and shear modulus (and thus wavenumber where ). The entire host domain of the composite is denoted by *D*_{0}.

Displacements in the *x*_{3} direction denoted by where are governed by the scalar wave equation(2.3)where is the indicator function(2.4)The Green function for the scalar wave equation in the host medium is governed by(2.5)and is given by(2.6)where is the (outgoing) Hankel function of the first kind and of the zeroth order. For ease of notation we drop the superscript (1) from all Hankel functions throughout. Thus, for orthotropic materials, we may restate the governing equation (2.3) together with the boundary conditions of continuity of displacement and normal stress on the interfaces between host and fibres in the form (Martin 2003)(2.7)where and are the deviations in mass density and shear modulus, respectively, and . Note that since the medium is assumed to be periodic and of infinite extent, there is no contribution from any external boundary in (2.7). The Green function is already non-dimensionalized and so on taking , where is a typical magnitude of displacement, and on dropping hats, the (non-dimensionalized) integral equation takes the form(2.8)in which and . In scaled coordinates, the lattice vector (2.1) becomes(2.9)where , and the cross-sectional area of the periodic cell is(2.10)

## 3. Integral equation methodology

We introduce an extension of the integral equation methodology initiated by Brazier-Smith (2002) and subsequently discussed and formalized by Parnell (2004). This scheme generates a hierarchy of equations, the first step of which is to differentiate the governing integral equation (2.8) with respect to () to give(3.1)where we have used the (anti)symmetry of derivatives of the Green function, , . Let us define the position vector of the origin of the fibre and suppose that . Upon taking the Taylor series expansion of the Green function about (as defined in (2.9)), we may write the governing equations (2.8) and (3.1) in the form(3.2)(3.3)where denotes the *i*th derivative with respect to . Note that the Green function to be summed can be treated as a continuous function of ** p** so that differentiation with respect to is well defined. Subsequently, it is summed over the discrete set of points defined by , which correspond to distinct lattice vectors . Furthermore, the terms involving require separate treatment owing to the singular nature of the Green function in this region. We have also introduced the

*displacement*and

*displacement-gradient moments*(of order ), defined as(3.4)(3.5)where the orders are deduced by noting that

*w*is a smooth function so that

*w*and all its derivatives will be . Recalling that is the position vector of , we use the notation to denote the operation of multiplying equation (*) by and integrating over the domain in the

**plane. Therefore, applying [(3.2)] and [(3.3)] for**

*x**k*=1,2 and again taking the Taylor series expansion of the Green function and its derivatives about

**=**

*x***(with the exception of the separated terms), we obtain(3.6)(3.7)where(3.8)Regarding the separated terms, it is clear that , and although inside , we find that because is singular in . Therefore, since and , the only separated term that contributes at , and thus at leading order, is(3.9)To solve this system (at leading order in**

*r**ϵ*), we see from the right-hand side of (3.6) that we must determine each , and for . We should therefore use the operator for every in order to derive an infinite system of equations for these unknowns. However, let us truncate this system to retain only , and and denote this as the system; it will also be known as

*the leading-order system in ϕ*for reasons that will become clear shortly. Furthermore, assuming that all fibres are identical, the system becomes(3.10)(3.11)In this we have defined the volume fraction (per unit span in

*x*

_{3}) , where is the fibre cross-sectional area and is the periodic cell area given in (2.10). Let us now define the variable where we remind the reader that

**and**

*p***are the th and th lattice points, respectively, and pose**

*r**effective plane wave solutions*of the form(3.12)in which is defined as the non-dimensionalized effective wavenumber (scaled on ) in the direction , the angle subtended from the axis. The system of equations (3.10) and (3.11) (for ) therefore becomes(3.13)(3.14)(3.15)where(3.16)and(3.17)We shall now discuss the evaluation of each of the individual terms in (3.13)–(3.15).

### (a) Calculation of *I*_{ij}: the sum to integral step

In order to compute each *I*_{ij}, we shall turn the double summation into a double integral. If the function to be summed possesses an apparent non-integrable singularity, then we must remove this singularity before this step is taken. Consider(3.18)where is the area of the periodic parallelogram (see (2.1), (2.2) and (2.10)) and assume that the function is integrable. We may write this as a double sum of integrals over the periodic parallelograms as follows:(3.19)On taking the Taylor series expansion of in each of the cells, about and , we are left with the result that regardless of the geometry of the lattice, *at leading order in ϵ*,(3.20)With the exception of , all contain apparent non-integrable singularities that must be removed before this sum to integral step can be implemented. Information regarding the periodicity of the lattice will be retained in the term that is removed (a simple lattice sum) and the remaining double summation can then be turned into an integral following the steps above. This leading-order integral term arising in in (3.17) is therefore *identical for all geometries*, for fixed **γ**.

First, let us calculate which from (3.17) is defined as(3.21)where . Since the function to be summed possesses only an *integrable* singularity at the origin, we may immediately employ (3.20) so that(3.22)where we have defined the polar coordinate system(3.23)Summation is implied over repeated subscripts in (3.22) and we note again that is the direction of propagation of the wave. We show in appendix A*a* that when , on introducing , (3.22) becomes(3.24)where, here and henceforth, a function with a scalar argument γ_{1} indicates the abbreviated notation of *I*_{ij} (γ(0), 0). In fact the result (3.24) is rotationally invariant and thus holds for all .

It turns out that we may determine the effective density and shear modulus via decoupled equations and so, as will be seen shortly, it is not necessary to calculate and . In order to determine and (defined in (3.17)), we must first remove their respective singularities at the origin. For notational convenience, we shall define and , and thus(3.25)On performing the sum to integral step by noting that, as ,(3.26)we may write these lattice summations as(3.27)where(3.28)and(3.29)where we have introduced the second form in terms of a new function for convenience, as we show later. Choosing , we show in appendix A*b* that(3.30)Therefore, on defining , from (3.27)–(3.29)(3.31)The remaining coefficient to be determined in (3.17) is which can be shown to be for any lattice structure by removing the singular term as above.

### (b) Treatment of the ‘singular’ term

The separated term defined in (3.16) and (3.9) possesses an apparent singularity at . In order to study this further, let us write(3.32)where is a circular domain2 of radius having its origin at , the location of the apparent singularity. Thus, we have separated into a singular part and a non-singular part on domains and , respectively. Since the singular region is excluded from , we may freely pass the derivative through this ** y** integral and interchange the order of integration in

**and**

*x***. On using the (anti)symmetry of differentiation of the Green function, we thus deduce that(3.33)where(3.34)The superscript**

*y**ϵ*in denotes its frequency dependence. We note that is directly related to , the component of the

*dynamic Eshelby tensor*associated with longitudinal shearing, by noting that (see Cheng & Batra (1999) for more details and Eshelby (1957), Mura (1982) and Michelitsch

*et al.*(2003) for more general commentaries)(3.35)This term incorporates the shape of the fibre cross section into our analysis and so we shall call the

*shape tensor*. Since , at leading order in

*ϵ*the shape tensor corresponds to static behaviour; thus, we can write(3.36)Furthermore, we can exploit the theory of the Eshelby tensor to state that will be

*uniform*(i.e. independent of

**) for inclusions of**

*y**elliptical*cross section (Mura 1982). Since in this article we restrict attention to the system, we need only retain the

*uniform*part of the static Eshelby tensor and thus we define(3.37)and note that for elliptical cross sections . Spatially dependent parts of in (3.33) will lead to higher order displacement-gradient moments that only appear in the full system,

*not*in the system. Therefore, on letting in (3.33) and noting the form of in (3.12), the contribution to the system for fibres of general cross section will be(3.38)For general cross sections, terms neglected are the extra contribution of higher order displacement-gradient moments (important only in the full system) and terms of higher order in

*ϵ*. For

*orthotropic*media, . The realization that this term is a scalar multiple of the Eshelby tensor is important since we may use the existing literature for inclusions of various shape (see Mura 1982). Similar terms arise in the full system, which correspond to weighted Eshelby tensors.

We now turn to the apparently singular part of , which has the form(3.39)Since is bounded inside this (small) region centred on , we may take a Taylor series expansion of the displacement gradients about in (3.39). This leads to integrals of the type(3.40)for and for . Referring to the definition of in (3.23), we may define the local polar coordinate system, centred on ** x** by and note that as (3.41)where are known constants. We may therefore write(3.42)It is straightforward to show that these integrals are zero either by virtue of the

*β*integral or by taking the limit . Therefore, we conclude that the apparently

*singular*term for any inclusion geometry.

### (c) The resulting {00} system for orthotropic media

On taking so that waves propagate in the *x*_{1} direction with wavenumber *γ*_{1}, the results of §3*a* and *b* infer that the {00} system is(3.43)in which(3.44)where(3.45)Note that we consider only orthotropic media so that . It is now clear from the form of (3.44) that *leading order in ϕ* means that the elements of ** A** contain only terms that are linear in

*ϕ*. The full system gives rise to a larger matrix with terms involving higher powers of

*ϕ*. For non-trivial solutions of (3.43), we set giving a condition on the effective wavenumber

*γ*

_{1}. From this we may deduce effective properties as will be discussed in §5. In order to determine the effective wavenumber in the

*x*

_{2}direction,

*γ*

_{2}, we could set and recalculate all the above but this is in fact unnecessary. What we can do instead is rotate the material itself by clockwise. This leaves the integrals evaluated as part of unchanged. We still consider waves propagating along the

*x*

_{1}axis, but we denote the wavenumber in this (new) direction as

*γ*

_{2}. The only difference is the resulting set of lattice summations and the position of the shape tensor in the system matrix. Thus,

*γ*

_{2}is determined by setting , where(3.46)Note in particular that if

*m*=1, and define (

*j*=1,2) as(3.47)Alternatively, if

*d*=1, it is defined by(3.48)

## 4. Lattice summations and shape tensors

The effective wavenumber has been obtained in §3 in terms of several basic quantities regarding the structure of the composite, including the shape tensor and the lattice summation. These must therefore be calculated for any given composite material and the objective of this section is to illustrate the solution procedure. This is done for circular–cylindrical inclusions on a rectangular or hexagonal lattice.

### (a) Rectangular lattices: orthotropy

Let us first restrict attention to a rectangular lattice with aspect ratio so that we can choose(4.1)and thus the periodic cell area is . As , for , the lattice sum (3.29) becomes(4.2)(4.3)which is rapidly convergent. On introducing the form (3.29), we see that for (4.4)where we note that the only difference between and is the position of the factor *A*. This comes about as a result of the order of summations in the double sum (4.2). In the case of , the wave-like term in (4.2) is . On the other hand, for , this wave-like term becomes due to the rotation of the material by . Anisotropy in the composite material is therefore induced by these lattice sums as well as the fibre cross-sectional shape. It turns out (provable by Keller's identity; see §5*b*) that these summations satisfy the relationship(4.5)and therefore in the special case of a square lattice (tetragonal symmetry) when *A*=1 we find that(4.6)

### (b) Hexagonal lattice: transverse isotropy

We now consider fibres positioned on a hexagonal lattice. The parallelogram lattice vector defining the positions of the fibres is given by(4.7)so that . It is easier however to evaluate the hexagonal lattice sum as the superposition of *two* rectangular lattice sums of the type in (4.1) with in each but with one offset from the origin by . We achieve this by writing(4.8)Proceeding analogously to the above, this reduces to(4.9)and thus we conclude that for the hexagonal lattice(4.10)Similarly, we find that , as should be expected from symmetry.

### (c) Shape tensor for the {00} system

Note that for the quasi-static effective properties, we need to determine only the component and thus can expand for small argument to find(4.11)Furthermore, since here we are interested in the system (leading order in *ϕ*), we need only the uniform part (i.e. independent of ** y**) of this, denoted , as discussed in the paragraph following equation (3.36) and noting equation (3.37). Therefore, by incorporating the shape tensor , a general result for the effective properties of the composite may be deduced for fibres of arbitrary cross section. Since we restrict attention to

*orthotropic*media, we consider only fibres with reflectional symmetries in the and axes. Therefore, we can write(4.12)Furthermore, from potential theory,(4.13)Since both and are scalars, we can write where so that, on performing a contraction on in (4.12) and exploiting (4.13), we see immediately that(4.14)Finally, we note that for fibre cross sections that have a rotational symmetry of (e.g. circles, squares) we must have and thus , so that(4.15)This latter point means that the result gives the same expression for fibres of square and circular cross section, the difference being of course that the

*volume fraction*will be different in each case. Furthermore (higher order) shape effects are taken into account when incorporating more equations in the truncated system, as described previously.

## 5. Effective material properties

For general directions , the form of for a homogenized orthotropic material having mass density (scaled on ) and longitudinal shear moduli and (with respect to the and axes, respectively, and scaled on ) would be(5.1)If we take *m*=1, then and we therefore find that(5.2)and similarly for *d*=1 we have and therefore(5.3)Thus, in order to determine the effective density we can simply take *m*=1, and to find the effective shear moduli we can take *d*=1 together with appropriate choices for . Alternatively, we showed above that we could first take in order to determine , then rotate the material by clockwise and consider for a second time to determine .

### (a) Effective density

On setting *m*=1, it is clear from (3.47) and (5.2) that the leading-order effective wavenumber is independent of and therefore defines the effective mass density of the composite as(5.4)which is the usual arithmetic average, being dependent only on volume fractions of the fibres and *not* on their shape.

### (b) Effective longitudinal shear modulus

On setting *d*=1, we see from (3.48) and (5.3) that for ,(5.5)Note the rich structure in the solution (5.5). Four separate terms account for different physical effects: the volume fraction ; the shear moduli ratio *m*; the lattice summations ; and the tensor incorporating the shape of the fibre cross section. We also note that this *leading-order system in* is not equivalent to a simple linear dependence on associated with dilute schemes. As we shall show shortly, the predictions of effective moduli given by this leading-order scheme will give excellent results even at relatively high volume fractions. Note that on employing (4.5), the definition of in (3.45) and the expressions (4.12) and (4.14), we see that in (5.5) satisfies the phase exchange identity (Keller 1964)(5.6)Alternatively, ensuring (5.6) holds for arbitrary cross section proves the result (4.5).

Restricting attention to circular cross sections, using (4.15) we find that(5.7)where . In the case of both square and hexagonal lattices, , thus giving the same result in terms of volume fraction of fibres. Clearly though for a given radius, the volume fraction of fibres in the two materials will be different. In figure 2 we plot the effective longitudinal shear modulus (scaled on ) for a square lattice material with and *m*=10, respectively. We show results for the integral equation method (IEM) discussed above, given by (5.7) with and also those given by the method of asymptotic homogenization (MAH) determined in Parnell & Abrahams (2006). With reference to this latter paper, it is in fact possible to show (via some rather lengthy algebra) that the results obtained in equations (3.25) and (3.26) of Parnell & Abrahams (2006), by classical asymptotic homogenization with only a *single term* in the multipole expansion described subsequently in equations (4.45) and (4.46) of that paper, are *identical* to the expression (5.7). We also plot energy bounds given by Hashin & Rosen (1964), which we note give a significant improvement on the Hashin–Shtrikman bounds by exploiting periodicity. Our main consideration is how close the leading-order integral equation scheme approximation lies to the fully converged asymptotic scheme. Results are plotted against the radius *r* of the fibre so that corresponds to touching fibres. We plot the corresponding results for a hexagonal lattice in figure 3. Finally in figure 4 we compare the predictions of the IEM with the MAH for a rectangular lattice with *A*=2 (so that the material becomes macroscopically anisotropic). Note that in all cases the leading-order integral equation scheme shows excellent agreement with the fully converged asymptotic scheme up to fairly high radii, thus providing a very simple and useful approximation of effective properties.

## 6. Conclusions

In this article we have introduced a new integral equation approach to homogenization and shown that its leading-order system gives predictions of static effective moduli for periodic fibre-reinforced media that agree extremely well with existing methods. The leading-order result in was given for fibres of arbitrary cross section and results were plotted for circular shapes.

We have indicated the rich structure of the solution for the effective longitudinal shear modulus where distinct terms account for specific physical features such as fibre cross-sectional shape, material properties, etc. This separation of physical phenomena is in contrast to the MAH (e.g. Parnell & Abrahams 2006) where all this detail is obscured by the fact that a doubly periodic cell problem must be solved; for fibres of arbitrary cross section, this requires some numerical scheme such as a collocation technique. This is not the case for the present integral equation formulation where all the details of the inner *local* problem and its periodicity are conveniently embedded in the shape tensor and the lattice summation terms. Both of these are easily determined; the former by recourse to the existing literature on weighted Eshelby tensors (via polynomial eigenstrains) (Mura 1982).

Higher order contributions can be found by incorporating more terms in the original linear system of equations involving displacement and displacement-gradient moments. This system was truncated at leading order in this article to obtain a result that retains the interaction between fibres (since ) and is therefore *not* equivalent to dilute schemes. The full system (which can be calculated algorithmically) will be discussed in a later article together with the results for fibres of more complex cross section. It turns out that the effective longitudinal shear modulus is modified from that given in (5.5) and takes the form (6.1)where the coefficients depend on higher order shape tensors (weighted Eshelby tensors), higher order lattice sums and the shear moduli contrast *m*. This dependence can be determined explicitly. We see from this form why the leading-order result (in ) obtained in this article for out of plane shear properties is so accurate. By retaining only the leading-order term, the error in the numerator and denominator in (6.1) is . Note that it may also be shown that retaining each additional order in in (6.1) is exactly equivalent (for circular cross sections) to the addition of an extra term in the multipole expansion in Parnell & Abrahams (2006), where expressions were obtained via the MAH. The form of (6.1) is, however, preferable for the reasons outlined above.

The method described has been termed an *elastodynamic* technique since we can determine higher order corrections in *ϵ* and thus predict dispersive effects when the frequency is increased (while retaining the notion of separation of scales). This is beyond the scope of the present article and results will be presented by the authors elsewhere.

It is quite clear that the method may be extended to fibre-reinforced media that possess more complicated macroscopic anisotropy such as monoclinic materials. Furthermore, phases that are themselves anisotropic could be discussed. This is another distinct advantage of an integral equation approach and again this will be discussed in the articles to follow. It would also be advantageous to develop the method for random media. Work is in progress on this problem and will be reported on in future articles. Furthermore, it is clear that the approach would be very useful for determining the effective properties of particulate composites in three dimensions. Although schemes have been developed successfully in this case by Iwakuma & Nemat-Nasser (1983), Nunan & Keller (1984) and Sangani & Lu (1987), for example, they are often inefficient and/or rather cumbersome. The rational approach of the present method, and the elegant and efficient form of its solution, offers an attractive scheme for complex composites.

## Acknowledgments

W.J.P. wishes to thank EPSRC and Thales Underwater Systems Ltd for their financial support while carrying out this work.

## Footnotes

↵Note that

*not all*lattice vectors of this form give rise to an orthotropic material, but all orthotropic materials may be described in this form.↵Clearly when

approaches the boundary of*x**D*_{ab}, this domain must be modified accordingly.- Received October 8, 2007.
- Accepted January 29, 2008.

- © 2008 The Royal Society