Variational and Related Methods for the Overall Properties of Composites

ADVANCES IN APPLIED MECHANICS. VOLUME

21

Variational and Related Methods for the Overall Properties of Composites J . R . WILLIS School of Mathematics University of Bath Bath. England

I . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I1. Preliminary Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A. Types of Composites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B. Definitions of Overall Properties . . . . . . . . . . . . . . . . . . . . . . . . . .

2

3 3 1 13

I11. Methods Related to the Classical Variational Principles . . . . . . . . . . . . . . . A. Small Variations in Moduli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B. Bounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV . Methods Related to the Hashin-Shtrikman Variational Principle . . . . . . . . .

13 20

A. The Integral Equation and Variational Principle . . . . . . . . . . . . . . . . . B. Perturbation Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C. Variational Estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23 25 33

23

V . Self-Consistent Estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

VI . Generalizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A. Viscoelasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B. Thermoelasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VII . Problems Which Lack Convergence . . . . . . . . . . . . . . . . . . . . . . . . . .

47

A. A Model Problem: Diffusion to a Random Array of Voids . . . . . . . . . . . B. An Integral Equation and Perturbation Theory . . . . . . . . . . . . . . . . . . C . Variational Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48 50 55 56 58 60

VIII . Wave Propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

A . Polarization Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B. Variational Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66 72

IX . Recent Developments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74 14

1 Copyright 0 1981 by Academic Press. Inc . All rights of reproduction in any form reserved. ISBN 0-12-002021-1

2

J. R. Willis

I. Introduction A wide range of phenomena that are observable macroscopically are governed by partial differential equations that are linear and self-adjoint. This article is concerned with such phenomena for materials, such as fiberreinforced composites or polycrystals, whose properties vary in a complicated fashion from point to point over a small, “microscopic” length scale, while they appear “on average” (that is, relative to the larger, macroscopic scale) to be uniform or perhaps, more generally, their properties appear to vary smoothly. The determination of these “overall”properties from the properties and geometrical arrangement of the constituent phases generally requires the solution of boundary value problems for randomly inhomogeneous media. The development of approximations is necessary for two reasons. First, even if the geometry of the medium were known exactly, the solution of any particular boundary value problem would be hopelessly complicated. Second, only certain features of the geometry (typically, concentrations of phases and low-order correlation functions) will be known in practice. In this situation, the best approximations are likely to be those that are derivable from a variational principle: this provides the unifying theme of this article. Problems are formulated in terms of operator equations that are self-adjoint and so generate associated variational principles. The equations are then solved, either by formal perturbation theory or by use of some ad hoc simplifying approximation, but usually the approximations are subsequently shown to have some precise variational status. Most of this article is concerned with static problems, for which many of the variational estimates of overall properties actually provide bounds. This article is concerned more with methods and principles than with results for particular media, though illustrative examples for which the algebra is not too complicated are given at a number of points. Whether or not a variational approach is adopted, a problem appears when any attempt is made to replace a finite body by one that is infinite: direct replacement of the Green’s function for the finite body by the infinite-body Green’s function leads to integrals that may converge conditionally in some problems and diverge in others and so cannot be allowed. The approach that is recommended is to start from a precise formulation for a finite body and to apply to this some mathematical analysis; this contrasts with an alternative approach, termed “renormalization,” which attempts to correct the offending terms without detailed reference to the finite-body problem. Much of this article treats the elastic behavior of composites, but it is emphasized that a number of other properties (conductivity, viscosity of a suspension, etc.) are described by the same equations. Extensions to visco-

Overall Properties of Composites

3

elastic and thermoelastic behavior are presented, for both of which the variational characterization given is believed to be new. Problems such as the resistance to flow of viscous fluid through a fixed bed of particles are mentioned and a model problem that involves diffusion is presented in some detail. This displays the same difficulty in relation to divergence of an integral and is the one problem of this type that has so far been approached variationally. In the concluding sections wave propagation is considered, using a formulation that is a natural extension of that used for static problems. This is under development at present and only approximate solutions of the operator equations have so far been obtained. However, this article is concluded with the derivation of a variational principle against which it should be possible to judge the approximate dynamical solutions. It is likely (but this cannot be guaranteed because it has not yet been done) that all of the procedures developed for static problems have their dynamic analogs so that, in the future, a unified dynamical treatment may be possible, which includes the static limit as a special case. 11. Preliminary Definitions

A. TYPESOF COMPOSITES For the purpose of this article, a composite is a particular type of inhomogeneous body, whose properties vary from point to point on a length scale (the microscale)which is small in comparison with overall dimensions. Generally, it will contain n distinct phases and perfect bonding will be assumed across phase interfaces; imperfections in bonding may, of course, be important in practice, but relate rather to the failure of a composite, a subject beyond the scope of this review. Examples of composites include glass fiber reinforced plastics, for which n = 2, the phases being the fibers and the resin, and polycrystals, for which n -+ co,each crystal orientation being associated with a different phase. It is reasonable for these two examples to treat the individual phases as continua and attention will be confined to materials (and phenomena) for which this approach is allowed. Thus, so far as the present work is concerned, a crystal such as MgO would be regarded as a homogeneous continuum. The problems that will be considered for composites may be described broadly as those in which the phases are not discernible individually, though they do, of course, determine overall properties. In the simplest case, the composite may appear macroscopically as a body whose properties are uniform, though more generally its overall properties may

4

J . R. Willis

vary: this would be the case, for example, if the mean concentration of reinforcing fibers varied over a specimen. The determination of these overall properties from those of the constituent phases is the subject of this article; they are defined more formally in the sections that follow. In preparation for this, a brief description is now given of the relevant geometrical properties of the microstructure. The simplest idealization of the microstructure of a composite would be to regard it as periodic, with period 6 much smaller than specimen dimension; for example, the composite may comprise a homogeneous matrix containing a regular array of inclusions. The properties of such a material would either be known with precision at each point or else might be subject to uncertainty only in relation to fixing the position of the inclusions relative to the boundary of the specimen: the uncertainty would be resolved if the position of one inclusion (and perhaps also the orientation of the lattice on which the inclusions are arranged) were given. Manufacturing processes are usually less obliging, however, and produce materials about which only statements “on average” can be made. It is desirable, therefore, to treat a composite as a random medium. For this purpose, any particular specimen is regarded as one with a label a which defines its configuration precisely. The parameter a belongs to a sample space 9’over which a probability density p(a) is defined. Then, any particular property f of the composite (such as the mass density at point x, or the product of mass densities at x and x’)is a function of a and has mean value, or ensemble average,

Of particular interest is the indicator function f,(x)that takes the value 1 if P,(x) of finding phase

x lies in phase r and zero otherwise. The probability r at a chosen point x is then expressible in the form

and the probability Prs(x,x’)of finding simultaneouslyphase r at x and phase s at x’is

Probabilities involving more points follow the same pattern. Mean values of other properties of the composite are expressible in terms of the essentially geometrical functions P,,P,,, etc. For example, if phase r has mass density p , and the density of the composite at x is denoted by p(x), then, in the case

Overall Properties of Composites

5

of n discrete phases,

Similarly, n

(P(x)P(x’)>=

n

C C PrPsPrs(X,X’)* r = l s=l

(2.6)

Correlation functions, of which (2.6) is an example, will be featured prominently in the analysis to follow. Reduced forms for them have been discussed by Miller (1969a,b) for a particular geometric model in which the space occupied by the composite is divided into cells, to each of which a value of r is assigned. However, the work contains more detail than is needed here, particularly for correlation functions of high order, and also restrictions (particularly isotropy) are imposed that will not always apply in the sequel. The form that probabilities such as Pr, take when phases are distributed isotropically has most recently been discussed by Fokin (1979). It should perhaps be mentioned that definitions like (2.2) to (2.5) rely only upon the notion of the random field f,(x) which need not be associated with a cell model, even though many composites may conveniently be viewed in that way. In the limiting case of a continuous rather than discrete variable r, the probabilities P r , Prs, etc. will be interpreted as densities and the sums in (2.4), (2.5) replaced by integrals. A more general approach based upon is a refinement that is measure (both here and for the sample space 9’) unlikely to be needed in practice; it is definitely not required in the present work. The foregoing discussion applies to any composite. However, for a composite which comprises a matrix containing a distribution of inclusions, there is some advantage in developing a description that allows for this feature. Elaborate notation is avoided by considering for illustration a matrix that contains a single set of inclusions, of identical properties, shape, and orientation. If the individual inclusions are labeled A = l, 2, . . . , the configuration of the composite is defined by the positions { x a ; A = 1,2,. . . } of any conveniently chosen point of each inclusion. For ease of description, the point x A will be referred to as the center of inclusion A, regardless of how it is chosen. In place of (2.2), it is now desirable to define a probability density P A for finding an inclusion centered at x A . This is achieved by evaluating the mean number of inclusions centered in any subset U of the region V

J . R. Willis

6

occupied by the composite, to give

where Y u ( k ) is the subset of Y for which there are exactly k inclusions centered in U . Similarly, the joint probability density PA, for finding distinct inclusions centered at x A , X , is defined so that, for disjoint subsets U and W of V,the expected value of the product of the number of inclusionscentered in U with the number centered in W is expressible as

where Y,,,(k, 1) is the subset of Y for which there are exactly k inclusions centered in U and 1in W .In the particular case where the composite is known to contain exactly N inclusions, Y v ( k ) = Y if k = N and is empty otherwise. Then, Eq. (2.7) implies Jv

dxA = N ,

(2.9)

as it should. Also, if U is so small that it can contain at most one inclusion, and if W = V - U , then Y,,,(k, 1) = Y,(k) if 1 = N - k and k = 0 or 1 and is empty otherwise. Equation (2.8) then implies

The conditional density PBIA is defined by the relation PAB= P B ~ A P A ,

(2.11)

with PBIA = 0 whenever inclusions A and B would overlap. It follows from (2.10) that (2.12) again as expected. The discussion so far covers random media more general than composites, in that no particular microscale has been assumed. Suppose now that the dimensions of the region I/ occupied by the composite are large compared with a typical microstructural dimension such as grain size or mean inclusion spacing. The notion of a statistically uniform material may be defined by the requirement that probabilities such as P,(x),Prs(x,x’), P A ,PA, should be insensitive to translations, so that P,, P A become constants, while P,,, PA, become functions of (x’ - x), (xB- xA),respectively. Of course, for any particular finite V , these requirements cannot be met exactly, but it is

7

Overall Properties of Composites

reasonable to assume their validity except in some “boundary layer” region close to the surface 8V of V . Another assumption that also requires V to be large will be made from time to time. This is that a composite possesses no long-range order. The mathematical implication for the joint probability Prs(x,x’) is that (2.13) when x and x‘ are far apart. This notion is independent of statistical homogeneity. An example of a composite for which (2.13) would not be realized is provided by a periodic structure, for which the only uncertainty may be in the precise location of phase r. Similar concepts may be applied to higher order probabilities and also to the densities P A , P A B ,etc., for a composite containing inclusions. Finally, for statistically uniform media, it is usual to make an ergodic assumption that local configurations occur over any one specimen with the frequency with which they occur over a single neighborhood in an ensemble of specimens. Subject to this assumption, the probability P , (which is now independent of x) may be estimated as the volume average of the function f,(x) for one specimen, P,, (which now depends upon x’ - x) is estimated as the volume average of f,(x x”)f,(x’ x”) (the integration being over x”), and so on. This assumption will be invoked at several points in the arguments that follow. It is as well to admit now, however, that precise conditions under which statistical uniformity and ergodicity(in the sense given)are approached asymptotically for large V are unknown at the present time.

+

+

B. DEFINITIONS OF OVERALL PROPERTIES The basic problem under consideration is to determine the response of a specimen or a structure when any kind of generalized loads (which could take the form of imposed forces, or displacements, surface temperature, electric field, etc.) are applied to it. To be specific, this initial discussion will be phrased in terms of the elastic stress analysis of a composite, since this particular aspect of composite behavior is considered more extensively than any other in the present work, but many of the remarks to follow will apply more generally. For the stress problem, then, for any particular specimen, stress components oijare related to strain components eijthrough the tensor of elastic moduli whose components are L i j k l , so that (2.14)

the strain components being related to displacement components uiby eij =

+ ujJ =

(2.15)

8

J . R. Willis

The summation convention for repeated suffixesapplies in (2.14). In (2.15) the suffixj denotes partial differentiation with respect to x j and the brackets applied to suffixes imply symmetrization. To avoid complicated formulas, a condensed notation will be used wherever possible so that, for example, (2.14) becomes u = Le.

(2.16)

In keeping with this shorthand, u, e, and L are referred to as stress tensor, strain tensor, and tensor of moduli, respectively. The tensor of moduli L varies over the region V occupied by the composite, taking the value L, in phase Y. If boundary conditions (tractions, displacements,or some mixture of these) are applied over dV and body forces fare specified over V ,the problem is to determine the stress, strain, and displacement throughout the composite. This is accomplished in principle by solving the equilibrium equation diva + f = 0

(2.17)

together with the constitutive relation (2.16) and the given boundary conditions. The moduli L vary over such a small scale, however, that the solution would provide an excess of detail even if it could be found and it would be desirable to extract from it some “locally averaged” information. This would consist, perhaps, of the values of stress, strain, and displacement, averaged over some representative volume, small on the scale of V , but large on the scale of the microstructure; such averages might be measured, in practice, with strain gauges or transducers. For some applications, it might also be appropriate to have some information on local fluctuations about the averages. For composites with periodic microstructure, these notions can be made precise by using the method of multiple scales (Cole, 1968), in which stress, strain, and displacement fields are expressed as functions of the macroscopic variable x/D and the microscopic variable x/6, where D, 6 represent length scales for the whole specimen and the microstructure, respectively. The fields are then represented as periodic functions of x/6, whose mean values and amplitudes vary with x/D. The approach was initiated by Sanchez-Palencia (1974) and forms the basis of the method of homogenization, comprehensively reviewed by Bensoussan et al. (1978). The method of multiple scales gives equations for the mean values that are consistent with an “overall” constitutive relation u

= Le,

(2.18)

where E is uniform and defined by a prescription that will be given later, together with equations defining the local fluctuation. Furthermore, the multiple scale solution has been proved to be asymptotic to the exact solution in the limit 6/D -+ 0.

Overall Properties of Composites

9

For the random composites with which this article is concerned, no corresponding precision has yet been achieved. There are, nevertheless, physically plausible prescriptions which are consistent with what is known for periodic composites. These are now outlined. In contrast to the periodic case, however, the prescriptions generate stochastic equations for which approximate solutions must be sought. 1. Dejnitions Based on Volume Averages

If it is assumed that a specimen behaves “on average” as though it were homogeneous, with constitutive relation (2.18), the overall moduli may be determined practically by a suitable set of tests in which stress and strain are uniform, in some nominal sense. To be specific, the boundary aT/ of the specimen could be subjected to the displacement u.1 = e..x. x E dV, (2.19) 1J 1’ which would generate the strain field 5 throughout V if the specimen were actually uniform. In fact, e varies over V , but its average over V may be expressed in the form

n

=

I/V J ; z i j d V = zij

(2.20)

by application of Gauss’ theorem, together with (2.19) on 8V. In (2.20) V is also used to represent the volume of the region V , to avoid introducing more symbols. Thus, when (2.19)is applied, the average strain over V is E precisely. Overall moduli E may now be defined by the relation

a = L5,

(2.21)

where denotes the average value of o over V ,when the boundary condition (2.19) is applied. It is interesting to note also that the average value 0 of the energy density over V may be calculated as (2.22) since aij,j= 0. Substitution of (2.19) into the last of these expressions and reversing the argument now gives (2.23) Thus, the moduli L may be defined either from the mean stress or from the mean energy density, if boundary conditions (2.19) are imposed.

10

J. R . Willis

Dually, the specimen could be subjected to boundary tractions compatible with a uniform stress field, so that -

aijnj = aijnj,

The equality 1/V Jv

(Tik

dx = 1/1/

(2.24)

x E aV.

Lv

aijxknj

dS

(2.25)

guarantees that cr indeed has average value 5 over V and overall moduli L may again be defined by (2.20) where, this time, the average value F of the strain is found from the solution of the boundary value problem. It follows also that (2.26) by substituting (2.24)into (2.22). Thus, again, overall moduli may be defined either from the average stress and strain or from the average energy density. These observations, made by Hill (1963a), provide rigorous definitions for overall moduli, for one particular specimen. They relate to tests that can actually be performed and apply whether or not the overall constitutive relation (2.18) is useful for predicting composite behavior when other boundary conditions are applied. Indeed, the two separate boundary value problems defined by (2.19)and (2.24)would yield different values for L, if applied to an arbitrary inhomogeneousmedium. It is plausible, rather than rigorously established, that Eq. (2.18) describes “average” composite behavior when arbitrary boundary conditions are applied and the extent to which the boundary conditions (2.19) and (2.24) generate estimates of L that are in agreement might be regarded as a partial check on the validity of the “overall modulus” concept. It may be remarked, finally, that these two definitions of L agree for the periodic composite in which, away from a boundary layer close to aV, they both generate periodic stress and strain fields with mean values 8, F, which Sanchez-Palencia(1974)used to define E over a single cell. The importance of Eq. (2.23) for the mean value of the energy is that it allows the construction of bounds for the components of L for any particular specimen. The principle of minimum energy states that, when displacements are prescribed over aV, the actual displacement field is that which minimizes the energy functional U* = (2V)-’ Jv e$LijkfeZdx

= (2V)-I

Jv e*Le* dx,

(2.27)

where e* is the strain associated with any displacement field u* over I/ that takes the prescribed boundary values. The simplest displacement u* compat-

Overall Properties of Composites

11

ible with (2.19)is

u* = z..x. V J

(2.28)

1

throughout V ;this leads to the inequality

0 I (2V)-'

s,,

ELEdx = (2V)-'E

sv

r= 1

L,f,(x)dxE,

(2.29)

where f,(x) is the indicator function defined in the preceding section. Thus, if the volume concentration of phase r is cr,Eq. (2.29) shows that Lv - L is positive semidefinite, where the estimate L,, for the overall moduli is the Voigt average (Voigt, 1889)

Lv =

c crLr. n

r=l

(2.30)

Dually, if traction boundary conditions are applied over aV, the complementary energy principle states that the actual stress field is that which minimizes the functional

U,* = (2V)-'

f,,o$!4ijkto~ld~ = (2V)-' J,,@*MU* dx,

(2.31)

where M is the tensor of compliances, inverse to L and @* is any selfequilibrated stress field that satisfies the given boundary conditions. The simplest choice is n* = 8, which yields the inequality U

where M is the overall compliance tensor, inverse to L and Mr is inverse to L,. Equation (2.32) shows that M, - is positive semidefinite, where

m

(2.33) or equivalently, L - LR is positive semidefinite, where LR is the Reuss average (Reuss, 1929) (2.34) The observation that the Voigt and Reuss estimates provide bounds for components of E is due to Hill (1952). Its importance lies not in the quality of the bounds themselves, but in the principle established by their existence: given the definition (2.21) of L, with the boundary condition (2.19), for

12

J. R. Willis

example, it follows without approximation that L, - L is positive semidefinite. The result is useful, however, only when some other information is employed in parallel. The definition of Lv involves only volume fractions and makes no allowance for the detailed distribution of phases. Therefore, at the very least, statistical uniformity must be assumed in addition, before a useful L is even defined. Also, a composite could be composed of isotropic phases, distributed anisotropically. In this case, Lv and LR would be isotropic, though L plainly would not be. Bounds for components of L are still obtained, however, provided that its overall symmetry is known from separate considerations. Implications of the Voigt and Reuss bounds for a composite containing parallel fibers have been discussed by Hill (1964). Better bounds, containing more statistical information on the composite, will be discussed in later sections. First, however, a different definition of overall properties will be given. 2. Dejinitions Based on Ensemble Averages Any given boundary value problem for a composite is defined by the equilibrium equation (2.17), together with the constitutive relation (2.14) and the boundary conditions. If the specimen to which the boundary conditions are applied is regarded as one sample from an ensemble, with the stress, strain, and displacement fields depend upon parameter c( E 9, CI and their means (a), (e), (u) could be sought. Then, at any chosen point, overall moduli L could be defined so that (a) = (Le) = €(e>.

(2.35)

For any given boundary conditions, there is no guarantee that the tensor f, does not vary over I/. Furthermore, if solutions at a chosen point x are examined, may be sensitive to the particular boundary conditions selected : an indication that this is so appears in Section III,A, where a formal perturbation solution generates f, as a nonlocal operator. If, however, it is believed that may be adequately approximated by a tensor of overall moduli, it may, as in the preceding section, be estimated by considering particular boundary value problems for which (a) and (e) are uniform over I/. It is, for example, plausible that ( e ) = e when the boundary condition (2.19) is applied, so long as the composite is statistically uniform. This depends, however, on the truth of the ergodic assumption described in Section II,A and will be discussed further when solutions are actually constructed. Finally, for such problems, a further application of the ergodic assumption leads to the ensemble averaged version of (2.26),that the mean

13

Overall Properties of Composites

energy density ( U ) at x is expressible in the form (U)

(2.36)

= +(ae) = +(a)(e),

which will allow the construction of bounds for f, by substituting into the classical energy principles suitable approximations to a and e. Kroner (1977) calls (2.36) the Hill condition. It should be remarked that the calculation of (a), (e), and (u) for any particular boundary value problem could be pursued as an end in itself, without explicit reference to f,. The drawback of the ensemble averaging approach is that mean values such as (a) have no strict meaning in any one sample, though it may be reasonable to regard (a) as some “locally averaged” version of a, provided that it varies slowly relative to the scale of the microstructure. The advantage of the approach, however, is that it can be followed even when (a) and (e) vary macroscopically and strict definitions based upon global volume averages cannot be applied. Usually, in the sequel, the ergodic assumption will be invoked and problems for which (a) and (e) are uniform will be emphasized. For wave propagation, considered in the final section, spatial variation is an intrinsic feature, however, and the ensemble averaging approach is the only one possible. 111. Methods Related to the Classical Variational Principles

A. SMALLVARIATIONS IN MODULI 1. General Perturbation Theory If the tensor of moduli L deviates only slightly from a constant value a formal solution to any given boundary value problem is easily developed in the form of a perturbation series. With the notation

Lo,

6L = L - Lo,

(3.1)

substitution of the constitutive relation (2.16) into the equilibrium equation (2.17)gives (LO)ijklUk,lj

+ (aLijklUk,J,j+ f;. = 0,

x

E

(3.2)

which is to be solved subject to a suitable boundary condition. Inessential complications are avoided by restricting the discussion to the displacement boundary condition u=uo,

x€av.

(3.3)

J. R. Willis

14

The problem defined by (3.2) and (3.3) can be expressed in the form of an integral equation by introducing the Green’s function G(x, x’) for a homogeneous “comparison” medium with moduli Lo. This satisfies (LO)ijklGkp,lj(X,x’) + 6ip 6(x - x’) = 0,

x E V,

(3.4)

Gip(x,x’)= 0, x E 8V. (3.5) If uo is now continued into V as the displacement field that the body force field f and the boundary condition (3.3) would generate in the comparison medium, it follows that the displacement u in the actual medium must satisfy (3.6) or, upon integrating by parts and using (3.5), Ui(X) = UP(X)

-J

aGip

7 (x, x’)8Lpjkl(X’)Uk,l(X’) dx’.

v axj

(3.7)

The symmetry of 6L ensures that only the strain e,, appears under the integral in (3.7). An integral equation for the strain field is produced by differentiation: symbolically, e = e0 - r&e,

(3.8)

where r denotes the linear operator

r:

z(x) +

Jv

r(x,x’)z(x’)dx’

(3.9)

whose kernel is related to G by rijpg(x, x’) = a2Gip/(dXjdXh)I(ij),(pq).

x’I-~

(3.10)

It should be noted that r has a singularity of order Ix at x = x’, which must be interpreted in the sense of generalized functions (Gel’fand and Shilov, 1964). It is convenient for present purposes to regard tensors of elastic moduli as operators with the form (3.9); with this interpretation, the kernel of 6L, regarded as an operator, would be 6L(x, x’) = 6L(x) 6(x - XI). Inversion of (3.8) gives e = [I

so that u = (Lo

+ r6L]-leo

+ &)[I + r6L]-’e0.

(3.1 1) (3.12) (3.13)

15

Overall Properties of Composites

The definition (2.32) off, relates (a) and (e); taking expectations of (3.12) and (3.13) and eliminating e0 shows that f, must be a nonlocal operator, given by

+ @L[I + r6~]-’)([1+

2 = L,

r6~1-l) -I.

(3.14)

The expectations in (3.14) are impossible to evaluate in closed form. However, when 6L is sufficiently small, formal expansion of the right side of (3.14) gives the perturbation series

2 = L~ + ( 6 ~ )+ (hL)r(&)

- (6Lr6L)

+ ’ . ..

(3.15)

For an n-phase material of the type described in Section II,B,l, the last term shown explicitly in (3.15) gives, when applied to (e), n

n

/

\

where 6L, denotes L, - Lo. Equation (3.16) demonstrates the nonlocal nature of the operator 2 and indicates, too, that successive terms in the series (3.15) contain integrals of joint probabilities for increasing numbers of points. Given these joint probabilities up to a certain order, an approximation to 2 follows by truncation of the series (3.15). Its evaluation is still not trivial because the finite-body Green’s function G from which is derived is generally not known. If V is large, it might be suspected that r could be replaced by r“,the corresponding operator derived from the infinite-body Green’s function G “ . This is not so, however, because G” is homogeneous of degree - 1 in (x - x’). Correspondingly, r“ is homogeneous of degree -3 and the integral in (3.16) is no more than conditionally convergent. One approach to the problem is to substitute I?“ into (3.16) and other terms in the series (3.15) and to single out for some special treatment those which involve integrals that do not converge. A “renormalization” is then effected, which amounts to identifying physically what the terms probably ought to contribute and replacing them by convergent terms that make this contribution. A discussion of renormalization is given by McCoy (1979). It will be apparent from the description just given that the author’s view of renormalization is that, while it might be necessary in contexts in which boundary conditions are unknown, it is no satisfactory substitute for straightforward analysis of a problem of which every aspect is fully specified.

16

J. R. Willis

An alternative derivation of is provided by the method of smoothing, given by Beran and McCoy (1970). In this, e is represented in the form e = (e)

+ e’

(3.17)

which, by substituting into (3.8) and averaging, produces (e) = e0 - r(GL)(e) - r(6Le’).

(3.18)

Subtracting (3.18) from (3.8) now gives the equation e’ = -r(6L - (6L))(e)

- T(6Le’

- (6Le’))

(3.19)

for e’, which may be solved by iteration to yield

e’ = {-I?(& - (6L)) + r p L r ( 6 L - ( 6 ~ ) ) - (6Lr(6L - (6L)))I - . . .>(e).

(3.20)

It follows from (3.17),(3.20),and (2.35) that f, = L,

+ ( 6 ~ )-

- ( 6 ~ ) ) )+ . . . .

(3.21)

This agrees with (3.151, but terms are automatically grouped so that r operates on a quantity whose mean is zero, wherever it appears. In analogy with (3.16), the last term shown explicitly in (3.21) gives, when applied to (e>,

c n

(hLr(6L

- ( 6 ~ ) ) )(e) =

n

r = l s=l

Jdxf 6Lrr(X,x ’ ) ( ~ , ~ xl) (x,

- P,(X)P,(X’))~L,(e(x’)),

(3.22)

in which the integral stands a reasonable chance of convergence if r is replaced by I‘“. McCoy (1979) regards (3.21) as a correctly renormalized series, but it must be pointed out that if, for a finite body, boundary conditions are applied for which (e(x’)) varies linearly (or worse) with x’, the integral converges only if P,, - PrP, approaches zero sufficiently rapidly at large separation of x and x’. The composite would thus at least be required to possess no long-range order: if p,, - P,P, merely oscillated about zero at large separations, as it would for a periodic composite, the integral would be conditionally convergent or even divergent for many choices of (e(x’)). When the probabilities have the required properties, however, Eq. (3.22) demonstrates that f, can be approximated by a local operator for any field (e(x’)) which displays negligible variation over the range within which P,,- P,P, is effectively different from zero. Furthermore, this local approximation may be obtained by considering a statistically uniform medium (for which P,, P,,,etc., are translation invariant and take values that agree with those given in (3.22)at the chosen point x), subjected to boundary

Overall Properties of Composites

17

conditions that generate a uniform field (e(x')). This restricted problem, for which a variety of other approaches are available, is taken up in sections that follow. The present section is concluded by remarking that the ideas summarized here have been developed and applied in many publications, amongst which those by Fokin and Shermergor (1969), Beran and McCoy (1970), Zeller and Dederichs (1973), Gubernatis and Krumhansl (1975), and Kroner (1977) are representative. Explicit recognition of the convergence difficulty is comparatively recent, and Zeller and Dederichs (1973) and Gubernatis and Krumhansl (1975) ignore it altogether. Fokin and Shermergor (1969) mention renormalization, but apparently mean by this the replacement of operators by formal perturbation series; no reference to convergence is made. Kroner (1977) combines r with a projection operator which has the effect of ensuring that it operates only on quantities with zero mean, as in (3.21). Earlier perturbation studies, for example those summarized in the book by Beran (1968), concentrate on solving the field equations (3.2) directly, though Karal and Keller (1964) performed firstorder perturbation theory, without considering convergence, on an operator equation like (3.8) for wave propagation. 2. Unform Mean Strain The general perturbation theory given above applies without change when the particular boundary condition (2.19) is applied. Then, if the composite is statistically uniform, it is to be expected that (e(x)> = c,

(3.23)

except, perhaps, in a boundary layer close to 8V and the local approximation to which f, reduces should be E,defined by (2.21) from the volume average iT of the stress. Some justification for these remarks is provided in the reasoning to follow, which derives essentially from work of Korringa (1973) and Willis (1977). Properties of the operator r are required and, as already anticipated in (3.9), it is convenient to set (3.24)

z = 6Le = (L - Lo)e.

The tensor z is called the stress polarization tensor. It will find extensive application later but, for now, it is simply a notational convenience. If v is defined so that (3.25)

v = u - uo,

it follows that v satisfies the field equations (LO)ijkluk,lj+ 7ij.j = 0,

E I/

(3.26)

J. R. Willis

18 with

v=o,

XEdV.

(3.27)

Equation (3.6) provides an integral representation for v, but with a view toward the ultimate relation of r to r", an alternative representation is now developed, which employs G" rather than G. It follows, in fact, by the usual argument involving Gauss' theorem, that ug(xr)=

fav -

G;(x - x')[tij

+ ( z i j ( x )- z$]njdS

jvG,POp,j(~- x r ) [ z i j ( x )- zi";.]dx,

(3.28)

where tij

=(WijklUk,l

(3.29)

and z* is any constant. Equation (3.28) represents v as a solution of (3.26) for any specified field z. The boundary condition (3.27) is also satisfied if T i = ( t i j + z i j - z t ) n j solves the boundary integral equation

Lv

G$(x - x ' ) T i ( x ) d S =

sv

GG,j(x - x ' ) [ z i j ( x )- z t ] dx.

(3.30)

Anticipating that the stress and strain n and e will oscillate rapidly about their mean values 8 and 8, z will likewise oscillate rapidly about its mean value 7 and the integrand on the right side of (3.30) oscillates about zero if z* is identified with 7.It follows directly from the field equations (3.26) that

L"

T,(x)dS = 0

(3.31)

and it is plausible that T i ( x ) actually oscillates about zero, when z* = f. Granting this, the boundary integral in (3.28) is significant only in a layer close to aV and so, elsewhere,

(3.32) Correspondingly, if the strain associated with v is denoted by e,

(3.33) the first equality coming directly from (3.8). Thus, when V is large in comparison with microstructure, Eq. (3.33) relates r to the infinite-body operator r"in a way that is consistent with (3.21). It should be noted, however, that the argument given is plausible rather than rigorous and would not apply at all for boundary conditions that would produce stress and strain fields that were not uniform "on average."

Overall Properties of Composites

19

Equation (3.33) approximates r by a translation-invariant operator that depends only on values 2, and so on the configuration of the composite, in a neighborhood of x', because the integral converges. The ergodic hypothesis therefore implies that, as x' ranges over V , the values of e(x') that are found occur with frequencies in proportion to the probabilities with which the associated local configurations occur. It follows that the average iZ of e over V may be identified with the ensemble mean (e) which may, in turn, be verified to be zero by taking the mean of (3.33). Therefore, (e) = V and now, from the mean of (3.12) with e0 = F,

F = ([I

+~ J L I - ~ ) ~ ,

(3.34)

so that this operator reduces to the identity when it is applied to a constant, though not otherwise. Equation (3.14) for f, now gives

fk- = h-, where

(3.35)

L = L~ + (GL[I

+r6~1-1) = L, + ( 6 ~ )- (6Lr 6 ~ +) . . . .

(3.36)

This result, although superficially simpler than (3.21), is, in fact, identical, since r is evaluated as in (3.33). This section is concluded by noting that perturbation theory may also be performed using elastic compliances rather than moduli. With the notation

M = Mo + 6M,

where

LoMo = (Lo + 6L)(Mo + 6M) = I,

(3.37) (3.38)

the integral equation (3.8) may be expressed in terms of u and uo by operating on each side with Lo + 6L to produce u = 60

+ (I - Lor)G L ( M+~ ~ M ) G .

(3.39)

AhMa,

(3.40)

It follows, using (3.38), that u =a '-

where

A

= L~ - L o n o .

(3.41)

The subtlety of the process of taking the limit of large V may be emphasized by noting that, as it stands, Eq. (3.40) applies when displacements are prescribed over V so that uo = LoE, which is different from a. It has been shown

J . R. Willls

20

already that (a) = a and so (3.40)gives -

L = ([I

+ A6M]-')Lo,

(3.42)

in which A is related to the Green's function for the displacement boundary value problem by (3.41).If, however, tractions were prescribed over dV, an equation like (3.40) could be derived, with A replaced by some operator A1 and 'a taking the value if. In this case, for consistency,

([I+A16M]-')a=if,

(3.43)

which shows that A, certainly differs from A. In fact, if A1 is used, Eq. (3.43) resembles (3.34) and it is plausible that, when V is large, A, takes the form

A1 : tl

+

LoFLo'l- Lo(tl- ii),

(3.44)

lr being defined by (3.33). If tractions compatible with a are prescribed over

dV, reasoning exactly parallel to that leading to (3.26) gives the expression

M=M,+(~M)-(~MA,~M)+-

(3.45)

for the overall tensor of compliances. It has, in fact, been verified that the first few terms of the series (3.42)and (3.45)agree, when A and A, are defined, respectively, by (3.41) and ( 3 4 9 , and r conforms to (3.33). Kroner (1977) introduced operators r and A for boundary value problems that were not precisely specified and related them through (3.41).Subsequent applications showed, however, that his A in fact was A l of the present work. Correct answers were generated because he employed the method of smoothing which has the effect of reducing A to A, in any case, as well as relating (a)and ( e ) independently of the exact choice of ao or eo. B. BOUNDS The formal expansions (3.36)and (3.45)for E and M are of limited practical use because they contain correlation functions of all orders. It was shown in Section 11, A how these correlation functions were related to joint probabilities for finding various phases at specified points, but only very few of these will be known in practice. Therefore, there will be no option but to truncate the series (3.36) and (3.45) after a few terms, to give conflicting estimates of L. If nothing is known about the composite other than the volume fractions c, of the phases, truncation of (3.36)at the term of first order in 6L gives the Voigt estimate -

L

'v

Lv = Lo

+ (6L)

(3.46)

Overall Properties of Composites

21

while a corresponding truncation of (3.45) gives the Reuss estimate -

L

-

LR = [Mo

+ (6M)]-’.

(3.47)

It was shown in Section II,B,l that the estimates Lv and LR at least have the precise status of bounding L from above and below and it will be shown now that other truncations of (3.36) and (3.45) similarly provide bounds. In fact, when e0 = e is prescribed, Eq. (3.8) has formal series solution cn

e=

C (-r 6 ~ ) ~ e .

k=O

(3.48)

Truncating after the mth term gives the approximation e* = B,F

(3.49)

where (3.50) The approximation (3.49) is associated with a displacement u* that satisfies the boundary condition (2.19) exactly. It is, therefore, a candidate for substitution into the minimum energy principle; this gives

ELCI (e*, (Lo + 6L)e*),

(3.51)

in which the “inner product” is defined by (3.52) Inspection of (3.10) shows that the operator r is self-adjoint so that, for any fields z1and t,,

(rz,,t,)= (r,,rz,);

(3.53)

this is also verified directly in Section IV,A. It also follows, from (3.4) and (3. lo), that

n o r= r.

(3.54)

The operators r,6L that act on the left member of the inner product in (3.51) may now be transferred to the right member, to yield

eLe I e(Lo + (6LBz,,,))e,

(3.55)

in which Lo appears by itself because Tz has mean value zero for any z. The ergodic hypothesis was also used to replace the spatial mean by the ensemble mean. The upper bound (3.55) contains correlation functions up to order

J . R. Willis

22

+

+

2m 1 and is precisely (3.36), truncated at 2m 1 terms. An exactly similar development, using (3.40) with goreplaced by 5 and A by A l shows, from the are obtained l by complementary energy principle, that upper bounds for % truncating(3.45)at an odd number of terms. These results are due to Dederichs and Zeller (1973). Improved bounds of odd order may be obtained by adopting, in place of (3.49), the trial field

e*

m

=

C (-r6L)kek,

k= 0

(3.56)

where ek are nonrandom “strain tensors,” arbitrary except that e, = e to satisfy the boundary condition. When m = 1, for example, substitution into (3.55) gives

~LizI s(Lo+ ( ~ L ) ) E- 2 ~ ( 6 L r&)el

+ e (6LT 6L + 6Lr 6LT 6L)e, ,

(3.57)

whose right side is minimized when

el

= (6Lr

6~

+ 6Lr 6Lr 6

~- l(6Lr ) b ~ ) ,

to generate the “optimized” third-order bound

€,

+ ( 6 ~ )- (6Lr 61,) x (ax6~ + 6Lr 6Lr 6

e,, where

(3.58)

= L~

~ - y) h L r 61,).

(3.59)

“Optimized” bounds of third order were first developed by Beran and Molyneaux (1966); the present account follows that of Kroner (1977). Bounds of even order have also been constructed by Kroner (1977) who has shown, essentially, that the last term in the series (3.55) is negative definite provided that Lois chosen so that 6L is positive definite: such bounds will be obtained from another standpoint in Section IV,C. Finally, it may be noted that the bounds given in this section improve as rn increases, since the trial form (3.56) includes the optimal form for any smaller value of m. Thus, the best bounds that can be obtained are those which use all of the correlation functions that are known. In practice,however, correlation functions or, equivalently,joint probabilities P,, . . .,may well not be available,though something may be known of their structure. In particular, isotropy of P , s , . can yield a major simplification in reducing the effect of the r operator to that of multiplication by a constant tensor P.Conditions under which this occurs for correlation functions up to a specified value of m define, in effect, Kroner’s (1977) concept of “disorder of grade m.” Kroner ,

23

Overall Properties of Composites used this in conjunction with trial fields of the form

[

e*= I - r &

1

n-1

C

e

k=l

(3.60)

to construct bounds for “disordered” materials. These will not be reviewed in detail, but bounds produced in Section IV,C for materials that need not be “disordered” reduce to some of those given by Kroner when disorder is also assumed.

IV. Methods Related to the Hashin-Shtrikman Variational Principle A. THEINTEGRAL EQUATION AND VARIATIONAL PRINCIPLE The integral equation (3.8) has so far been regarded as an equation for the strain e, but, as anticipated in the definition (3.24) of the stress polarization t,it can also be regarded as an equation in this variable; explicitly, so long as L - Lo is everywhere nonsingular,

[(L - Lo)-’

+ r]t = eo.

Small variations in moduli are not necessarily envisaged and the tensor of moduli Lo relates to any conveniently chosen homogeneous comparison material. Before considering the solution of (4.1), some of its properties will be noted. First, as indicated in the preceding section, the operator r is selfadjoint. This can be verified directly by employing the “virtual work” equality

ez) = 0 (44 for any stress cl that is divergence-free and any strain e2 derived from a displacement that is zero over i3V. Now let (01,

el = -rt1,

e, =

-rt2

(4.3)

for any chosen tl,t 2 ,and let b l be the stress associated with e l , so that cl = Loel

+ zl.

(4.4)

Then, ( z , , W = (Loel - bl, e2) = (Loel, e2), (4.5) by (4.2). The last expression in (4.5) is symmetric in the suffixes 42 by the symmetry of the elastic moduli Lo ; hence, r is self-adjoint. Next, setting

J . R. Willis

24

r is positive definite, and hence that (L - Lo)- + r is positive definite, so long as (L - Lo)-' is positive definite at each point

z1 = z2 in (4.5) shows that

of I/. If Lo is so restricted, therefore, Eq. (4.1) is equivalent to the minimum principle 2( 0, - 0 )I (z*, (L - Lo)- 'z*)

+ (z*, rz*)- 2(z*,eo)

(4.6)

for any approximation z* to z. The left side of (4.6) represents the extreme value of the functional -(z, eo) = (L,e

- a,eo) = 2(B, - 0),

(4.7)

where 2 0 , = (L,e, eo) = (e, L,eo) = (ao,eo) 2 0 = (a,eo) = (a,e),

(4.8)

(4.9)

by repeated application of (4.2). The functional on the right side of (4.6) is extremized by the solution of (4.1), regardless of whether L - Lo is positive definite. The variational principle is, in fact, a maximum principle when L - Lo is everywhere negative definite. To prove this, it is convenient to introduce a "strain polarization" 11 by z = Lev, 11 = MOT, (4.10) so that (4.11)

11 = M,a - e.

It now follows, using (4.2) where necessary, that = (11, Loti)

(z,

- (0,Moa).

(4.12)

But also, (z, (L - Lo)-lz) = (Lo%(L - Lo)-lLoll)

(L - Lo)- 'Lo

= Mo(M0 - M)-

- I.

(4.13) (4.14)

Hence, (z, (L - Lo)-lz) =

(119

(M,

- W'11) - (11,LoV)

(4.15)

and, by addition of (4.12) and (4.15), (z, (L

+

- ~ , ) - l t ) (z,rz) = (11, (M, - ~ 1 - 1 1 1 ) - ( a , ~ , a ) . (4.16)

This is negative definite so long as Mo - M or, equivalently, L - Lo, is negative definite; in that case, the minimum principle (4.6)is replaced by the

Overall Properties of Composites

25

maximum principle

2(u0- rli) 2 (t*,(L - Lo)-%*) + (z*,Tt*)- 2(t*,e0).

(4.17)

These derivations were given by Willis (1977),following arguments presented by Hill (1963b). Hill showed, in fact, that the Hashin-Shtrikman principle which is embodied in (4.6) and (4.17)can be derived from the classical energy principles, by discarding a term known to be definite: the even-order bounds of Kroner (1977)referred to in the preceding section are obtained by the same method and will later be derived from (4.6) and (4.17). The principle was originally derived, directly from the field equations, by Hashin and Shtrikman (1962a). B. PERTURBATION THEORY The case of small variations in moduli has, in effect, been adequately covered in Section III,A. It is useful to record, however, that perturbation theory for e, together with the definition t = 6Le, yields the formal series solution t = 6L

m

C (-r&)keo.

(4.18)

k=O

The stress polarization t plays a more distinctive role if a composite is considered which comprises a matrix containing a distribution of inclusions. For simplicity, the inclusions will be taken to be identical in shape, size and orientation and will be assumed to have moduli L1.The moduli of the matrix will be taken as L2 and the inclusions will be distributed at number density n, and associated volume concentration cl. Generalization to the case of several types of inclusion presents no more difficulty than the need for more cumbersome notation. If the comparison material is now identified with the matrix, so that Lo = L2, the definition (3.24) shows that t is nonzero only over the inclusions. Equation (3.8), with t = 6Le, defines e in terms of z anywhere in V , but the integral equation (4.1) now applies only over the regions occupied by the inclusions. It is convenient to call these regions V,, A = 1,2,, . . , as in Section II,A and to define't to be the restriction of z to V,. Equation (4.1) then gives [(L,

-~

~

)

+ r)tA+ C rzB= eo, -

l

B+A

xE

v,,

A = 1,2,. . . . (4.19)

The relation (0)

= Lo(e)

+ (.>

(4.20)

26

J. R . Willis

encourages study of the expectation value of z. This is expressible in the form (4.21)

where zi(x) denotes the expectation value of z,(x), conditional upon an inclusion being centered at x,, and V(x) is the set V(x) = (x,:

(4.22)

x E V,}.

An equation for z<(x) may be obtained by taking the expectation of (4.19), conditional upon V, being fixed. This gives

where '7: denotes the expectation of ,'t conditional upon inclusions being centered at xA and x'. Equation (4.23) thus introduces a new unknown, for which an equation could be generated by taking expectations of (4.19) with V, and V' fixed. This new equation, however, would contain &,, the expectation of zC given V,, V,, and V,, and continuing the process would generate an infinite hierarchy of equations for expectations of this type. The difficulty is resolved if a series solution is sought in powers of the concentration cl, in the limit c1 << 1. The integral in (4.23) may be expected to be of order c1 and its neglect generates the simple zeroth-order equation

[(L, - L ~ ) - I + r]zi

= eO,

x E V,,

(4.24)

which describes the response of a single inclusion V, perturbing a field that, in its absence, would take the value eo. This problem is of fundamental importance and is discussed below. 1. T h e Single Inclusion Problem

If the inclusion A is close to the boundary aV of V , use of the finite-body form of r is essential and little can be deduced from (4.24). If, however, the inclusion is far from aV, since it is isolated, the difficulties mentioned in Section II1,A for the operator disappear, and r may be replaced directly by r",whose kernel is given by (3.10)with G replaced by G". The infinitebody Green's function G" satisfies (3.4) with Lo = Lz and the boundary condition (3.5) is relaxed. G" is a homogeneous function of degree - 1 in x - x' and may determined by placing the point load at the origin, so that x' = 0. A useful representation is obtained by first noting the plane-wave decomposition (4.25)

27

Overall Properties of Composites

of the delta function, given by Gel'fand and Shilov (1964). This motivates study of the equation

+

sl,(g

(LZ)ijkfGfp,fj

*

x) = 0,

(4.26)

XI,

(4.27)

in which GCmay be taken to depend only upon (c * x). A particular solution is given by

G6 = - [LAC)] -

*

where L,(<) is a matrix with components (4.28)

= (LZ)ijkfljcf.

[LZ(C)lik

The superposition indicated by (4.25), applied to Gr, now gives G" in the form

Jla=

G"(x) = (8n2)-'

The integral equation (4.24), with inclusion is considered, becomes (Ll

- LZ)
7:

- (8n2)-1

jvAdx' 8°C

1

dS[Lz({)]-'S(C

(4.29)

x).

now replaced by z since only one

it,=

dSILZ(<)]ik-lcjlf

(X - x')]Tkl(x')l(i,j)= f?:(X),

XE

v,.

(4.30)

Equation (4.30) has a remarkably simple solution for the important special case in which V, is an ellipsoid and e0 is uniform over V,. Suppose first that V, is a sphere with radius a, centered at the origin. If z is taken constant over V,, the integral with respect to x' in (4.30) requires the evaluation of az/apZ

<-

J

VA

dx'qp -

c - xi),

when p = x. The integral over V, gives the area of the disc formed from the intersection of the plane * x = p with V,: this is n(a2 - p z ) if IpI < a and zero otherwise. When x E V', therefore, its second derivative is - 2n and the left side of (4.30) is constant over V,. Thus, when e0 is constant over V,, t is constant and given by [(L, - L 2 ) - l + P]t = eO, (4.31) where the constant tensor P has components 'ijkf

= (4n)-

it,, =

dS[L2(<)1,

'
(4.32)

A similar result is available when V, is the ellipsoid xTATAx< a. Evaluation of the integral is reduced to the case just considered by the dual transformations y = AX, 5 = ATC (4.33)

J . R. Willis

28 and produces the result f‘ijkl

J

= (4~lAl)-l

151=1

dS[LA<)li ‘tjtl[tT(ATA)-1 < ] - 3 ’ z l ~ i j ~ ~ k ~ ~ .(4.34)

Eshelby (1957) gave an explicit solution for the ellipsoidal inclusion problem when the matrix is isotropic, using methods of potential theory. He also proved that the field within an ellipsoidal inclusion in an anisotropic matrix would be uniform, but gave no detailed prescription for its determination. The earliest published derivation of (4.34) appears to be that of Khatchaturyan (1967), though (4.32) was given by Kneer (1965). Eshelby (1961) showed that, for an inclusion in an isotropic matrix, a polynomial field eo would generate a polynomial field within the inclusion. This was proved for an anisotropic matrix, in an unpublished essay by Willis (1970) and independently by Asaro and Barnett (1975);an application of the result was presented by Willis (1975). In general, evaluation of the tensor P requires a computation. Although progress can be made analytically for a number of special cases, the only one of sufficiently general interest and simplicity is that of a spherical inclusion in an isotropic matrix. If L2 is isotropic, therefore, with Lam6 moduli A, p, (L2)ijkl

= lzsij

6kl

+ p ( 6 i k sjl +

(4.35)

sjk)

and it follows without difficulty that (4.36) For a spherical inclusion, P is given by (4.32);it is isotropic and so may be represented in terms of a “bulk modulus” xP and “shear modulus” pp, so that Pijkl

= K p 6ij

+ pp(bik + 6jl 6jI

6jk

- S b i j 6kl).

(4.37)

The “moduli“ xP,pp may be obtained from the scalars f‘jikk

= 9Kp,

Pijij = 3 K p

+ 1oPp3

(4.38)

for each of which the integrand of (4.32) reduces to a constant. The result is expressible in the form 1

(3KP9

2PP) = (,(A

where K

=1

+ 2p)’ 15y::p))

=

(A. + )

3(K + 2pL) ’ (4.39) 4p)

5p(3k.

+ 2pj3 is the bulk modulus of the matrix. The notation p = ( 3 K P , 2PP)

(4.40)

29

Overall Properties of Composites

was introduced by Hill (1965a).It has the property that, for isotropic tensors A and B, their product AB may be represented as

AB = ((3K.4)(3KB), (2pA)(2pB)).

(4.41)

The identity has the representation (1, l), so that

A - ' = (1/31tA, 1/2pA).

(4.42)

If the inclusion is also isotropic, say with L1 = (3~',2p'),Eq. (4.31) has solution z = TeO,

where 3 K , = ([3(K'

- .,I-'

(4.43)

+ 3Kp}-1,

(4.44) 2.k = P ( P ' - P r l + 2 P P V This applies, of course, only when e0 is constant over the inclusion. Reverting now to the composite problem, overall moduli L may be estimated by choosing e0 = V. Then, ( e ) = C! and E may be deduced directly from (a), as given by Eq. (4.20). To lowest order in cl, (z) = c l z j since, when z i is constant, the integration in (4.21)simply multiplies the integrand by the volume of one inclusion. Thus, to first order in cl,

L = Lz + c,T.

(4.45)

2. Interactions In considering interactions between inclusions, it is helpful to rewrite Eq. (4.23)in the form C(L,- L J - ~+ rl.2 -tJrlx,(~,I, - p,)rz;

+ Jrlx,p,,,r(t.;B

-z ~ ) (4.46)

whose right side is precisely ( e ) , or v if e0 takes this constant value. Estimation of L to order c: requires an estimate of z$ to order c l r which, in turn, requires 0(1) estimates of z i , z;~ for substitution into the integrals in (4.46). z i is known already to lowest order, from the preceding section, while the zeroth-order estimate for zBABis defined by the equation (L, - LZ)-'zjB+ r z j ,

+

= eO,

x E V,

(4.47)

and a similar equation with A and B interchanged; these are obtained by taking expectations of (4.19)with V, and VBfixed and neglecting the integral, ~ ~(4.47) . describes the problem of two inclusions which involves ~ 5Equation

30

J. R. Willis

perturbing a strain field eo. For spheres in an infinite isotropic matrix and e0 = V, it was solved approximately by Willis and Acton (1976). More recently, Chen and Acrivos (1978a,b) have solved the same problem from a different standpoint. Both pairs of authors have produced estimates for the coefficient of c: in the expansion of L. The integrals that appear in (4.46) converge when r is replaced by r";Willis and Acton (1976), and subsequently Willis (1980a), whose formulation is otherwise followed in the present account, employed the explicit representation (3.33) for r from the outset. Chen and Acrivos (1978b) found it necessary to renormalize a conditionally convergent integral, which they attempted by a method successfully developed by Batchelor and Green (1972) for the related problem (formally the incompressible limit of the present one) of the overall viscosity of a dilute suspension. They noted, however, an ambiguity in the renormalizing prescription, which allowed at least two different estimates for L; an explanation of why a particular one (namely, the one consistent with the present development) is correct, formed a major part of their effort. The exact solution of (4.47) will not be considered, but before leaving this problem, a simplification of (4.46) will be discussed because it is relevant later. When the inclusions are spheres of radius a and PBIA is isotropic, and so a function of - xAl only,

,xI

(4.48) for any constant T;, where the tensor P is defined by (4.32). This may be seen by choosing, without loss, xA= 0 and expanding the left side of (4.48) explicitly in the form

(4.49)

XI

(,xI

< 2a is P. The integral over Now lx" - < 2a and the integral over > 2a is zero because PBIA - P, depends upon only and the integral of T"(x, + x" - x) over the surface = r is zero for any r > lx" this follows by differentiating with respect to r the corresponding volume integral over lxBl < r. A corresponding result applies for ellipsoidal inclusions, provided that PElAnow has "ellipsoidal" symmetry so that the composite is, in effect, generated from an affine transformation of an isotropically distributed suspension of spheres; an explicit proof is given by Willis (1978).

Ix,~

Ix,~

Ix,~

XI;

31

Overall Properties of Composites

Finally, Eq. (4.21) shows that ( 7 ) depends only on the mean value 5; of 7; over inclusion A. Integrating Eq. (4.46) with respect to x over V ’ reduces the kernel of the J? operating on 7; to the constant P which, in turn, generates P times the integral of t; over V,. Equation (4.46) therefore implies, to first order in el,

+ PI$

[(L, - LJ1

+ ClP?; - J ~ x B P B ~ S , ( &

=E

- 7;)

(4.50)

where, on the right side of (4.50), 7; and zBAB are required only to zeroth order and represents the mean value of r as x ranges over VA. 3. Closure Assumptions

At high concentrations el, the power series approach outlined above is not helpful. An approximation, whose status will be further clarified later, is obtained by assuming, in (4.23) or (4.46), that 7 BA B

- z BB

(4.51)

exactly. It is related to perturbation theory in the sense that, if the composite has no long-range order, equation (4.51) is asymptotically true at large separations, and few inclusions are close together when c1 << 1. Also, if the composite is perfectly ordered, with the inclusions arranged on a lattice, Eq. (4.51) is true identically at any el. Accordingly, (4.51) is called the quasicrystalline approximation. It was introduced, in the context of scattering theory, by Lax (1952). When it is adopted, Eq. (4.46), with e0 = F, reduces to

[(L, - LZ)-’

+ r]tj + J d x B ( p B 1 A

-pB)rt;

= s,

x E &, (4.52)

which retains some allowance for interactions through the integral. Equation (4.52) has a simple exact solution if the inclusions are spheres, distributed isotropically, or else are ellipsoids, with corresponding “ellipsoidal” symmetry for P B I AIt. is consistent to take 7; constant over VBsince the result (4.48) shows that the integral is then constant over V, . Therefore, (4.52) has solution 7;

= [(L,

-L

p

+ (1 - c,)P]-’e

and, correspondingly, L is estimated as

€ = Lz + C,[(L1

- LJ-1

+ (1 - c,)P]-’E.

(4.53) (4.54)

It may be noted that (4.54) agrees exactly with (4.45) to first order, but differs from the exact series for L at order c:, the discrepancy corresponding exactly to the integral involving & - 7; in (4.50). Although the limit c1 = 1

J. R . Willis

32

is not allowed, it is interesting that substituting el = 1 into (4.54) yields L = L1, so that there is at least some hope that (4.54) provides a reasonable estimate of L at high concentrations. The other important feature of (4.54) is that it contains no explicit dependence on the conditional probability though its isotropy (or, more generally, “ellipsoidal” symmetry) appears implicitly, through the tensor P which depends upon the matrix A defining the basic ellipsoid. The approximation (4.51) is one of a class whose next simplest member is

-

(4.55)

this provides a good approximation, except when three inclusions are close together and could be used to close the hierarchy of equations at the next stage, giving a closed integral equation for ~ 5 This, ~ however, . is not considered further. Closure assumptions like (4.51) or (4.55) may also be applied to a general n-phase composite. If the restriction of z to the rth phase is denoted by zr, n

(4.56) (4.57)

where z:(x) now denotes the expectation value of z’, conditional upon phase r being present at x. Substituting (4.53) into (4.1) and taking expectations conditional upon phase r being present at x gives (L, - ~ ~ ) - 1 z+;

f: J d x ‘ q x , x~)z;s(x~)~slr(x~, x) eo, =

s= 1

r = 1,2,. . . ,n

(4.58)

where zss(x’)denotes the expectation value of zs at x’,given the presence of phase r at x and phase s at x‘and PSlr(x’, x)is the probability for finding phase s at x’,given phase r at x. The analog of the quasicrystalline approximation (4.51) is zSs(x’)= z:(x’),

(4.59)

which reduces Eqs. (4.58) to a closed set for the unknowns t:(x). If e0 is uniform, and equal to F, the operator I? takes the form (3.33), except in a boundary layer close to LJV and it is consistent to take z: constant if the medium is statistically uniform. In that case, the equations reduce to the

33

Overall Properties of Composites algebraic set

in which the probability P, has been identified with the volume concentration c,. The integral is independent of x because PSI,is insensitive to translations. As for the suspension of spheres, further progress is possible if all of the functions PSI,are isotropic. The observation that the integral of r m ( x - x') over any spherical surface Ix = r is zero allows the integral in (4.60)to be replaced by one over an arbitrarily small sphere centered at x, over which PSI,- c, may be replaced by its value at x. Now

x'I

(4.61)

P,lr(X, x) = 6,s

and so, since the integral of Eqs. (4.60)simplify to

[(L, - Lo)-'

rm over the volume of

+ PIT:- P(T) =

e,

the small sphere is P,

r = 1 , 2 , . . . , n.

(4.62)

It follows from (4.62)that T: = T,(C

+ P(T))

where T, = [(L,

- Lo)-'

+ PI-'

(4.63) (4.64)

and hence that (2)

where (T)

=

= [I - (T)P]-'(T)E,

(4.65)

c:= c,T,. An estimate Z for L is then obtained in the form Z = Lo + [I - (T)P]-'(T).

(4.66)

If the material is in fact composed wholly of phase r, so that c, = 1 and c, = 0 for s # I, it is easy to verify that (4.66) gives f, = L,, as it should. C. VARIATIONAL ESTIMATES

It was shown in Section IV,A that the displacement boundary value problem is solved by the stress polarization z that extremizes the functional

+

F(T)= (T,(L - Lo)- lz) (z,r T ) - 2(t, e').

(4.67)

34

J. R . Willis

Furthermore, its extreme value may be expressed, when e0 = 5, in the form -

2(D0 - D)= Z(Lo - L)Z,

(4.68)

by (2.23). Thus, a variational estimate for L is obtained by substituting any trial field z* into (4.67). The first trial field that will be considered is Z* = 6~

rn

1 (-r6L)ke,

k=O

(4.69)

which is just (4.18), truncated at the mth term. The resulting F(r*) simplifies by use of the self-adjoint property (3.53) of r, and 6L correspondingly, to the form F(T*)=

- (Z,6LBzrn-I?),

(4.70)

where B, is defined by (3.50). The approximation (4.70) to F(T) thus leads, upon replacement of the spatial mean by an ensemble mean, to the estimate

Z = Lo + (6LBz,-l)

(4.71)

for L, which involves correlation functions for up to 2m points and is precisely (3.36), truncated at this order. The results of Section IV,A show, in addition, that - L is positive definite if Lo is chosen so that L - Lo is negative definite at each point of V. In this case, the estimate (4.71) becomes one of Kroner's (1977)bounds of even order. Dually, (4.71)yields a lower bound for L if Lo is chosen is that L - Lo is positive definite; this result was not given by Kroner (1977). He gave, instead, an upper bound of even order for the compliance tensor R. This could be obtained (together with a new lower bound) by working with the strain polarization tf defined by (4.10) and an integral equation like (4.1) for the traction boundary value problem, with I' replaced by A l : the Hashin-Shtrikman principle was given explicitly in terms of t f , as well as z, by Hill (1963b). 1. Hashin-Shtrikman Bounds

The next obvious development is to construct "optimized" bounds by substituting into (4.67) approximations z* that contain parameters. Here, the Hashin-Shtrikman principle displays a clear advantage over the classical principles discussed in Section III,B for, whereas the classical principles require either compatible strain fields or self-equilibrated stress fields, the polarization T* is subject to no such constraint. Considering first an n-phase material, it is natural to choose r* to have the piecewise-constant form (4.72)

35

Overall Properties of Composites

where the constants z, are arbitrary. Equation (3.67) then gives, when e0 = Z,

F(z*) = v-’

zr(Lr- Lo)-’zr Jdxf,(x)

r=l

- 2v-’

r= 1

(4.73)

zr Jdxf,(x)E.

Now by definition, the mean value of f , ( x ) over V is cr, the concentration of phase r. Also, except when x is in a boundary layer close to aV, the operator r takes the translation-invariant form (3.33) and so, when I/ is large,

v- 1 J d x ~ . ( xJdxrr(x, ) x~)fs(xf)

J d x / ’ r y x y f s ( x+ x y The mean value over x of f,(x)f,(x+ x”) is Prs(x,x + either by definition V - 1 Sd.r,(x)

CJ.

XI’),

of an “empirical” two-point probability or by use of the ergodic hypothesis directly. Hence, when V is large,

qz*) =

C,~,(L, r=O

- L , ) - I ~ ,+

cc n

n

r = l s=l

+

x (Prs(x,x x”) - crcs)zs- 2

t, J d x ” r y x ~ ~ ) n

1 crz$,

r= 1

(4.74)

which is independent of x because P,, is assumed to be insensitive to translations. This derivation has been given in detail as an explicit example of the ergodic assumption. Earlier, in Section III,A2, reference was made to the sampling by r of local configurations, but precise specification [such as (4.72)] was not possible and the argument necessarily remained qualitative. The best estimate for F(t) is now obtained by extremizing (4.74). This requires that (L, - Lo)-%, +

i:

s= 1

~ d x ” r ” ( x ~ ~ ) ( P , , r (-x c,~)z,~ , O=)e,

r = l , 2 , . . . , n,

(4.75)

having used the relation P,, = Psircrand chosen x = 0. It is interesting to note that Eqs. (4.75) are identical to (4.60), with z: = z,. Furthermore, the extreme value of F(z*) is simply - C:= crz$ and Eqs. (4.60), which were derived on the basis of the ad hoc closure assumption (4.59), in fact have the

36

J. R. Willis

status of providing the “best possible” estimate for L, relative to the HashinShtrikman variational principle. Finally, if this estimate is called L is positive or negative definite, whenever Lo is such that L, - Lo is correspondingly definite for all r. The bounds were given, at this level of generality, by Willis (1977).The original Hashin-Strikman bounds were derived rather differently from the variational principle by Hashin and Shtrikman (1962a,b, 1963) and, by a different method again, by Walpole (1966a,b), for the special case in which PSI,is isotropic and so a function of Ixrr/only. The integral operator in (4.75) then reduces to the constant tensor P, as in (4.62), and is given by (4.66). Willis (1977) noted, in fact, that (4.62) and (4.66) apply to a composite with “ellipsoidal” symmetry, if P is given by (4.34), which allows, for example, the construction of bounds for bodies containing aligned needles of finite length, or parallel cracks, as limiting cases. Kroner (1977) produced a form of (4.66), with the inverse operator expanded as a series, by substituting into his even-order bound the approximating field (3.60), taking the limit m -+ cc and simplifying the result by appeal to the concept of “disorder”; this appears to be realized as a limiting form of the cell model of Miller (1969a,b). The work of Hori and Yonezawa (1974) contains a relevant discussion. In view of the fact that “disorder” implies at least isotropic statistics and produces Hashin-Shtrikan bounds of less general application than (4.66), it receives no detailed treatment here.

z,

2. A Matrix Containing Inclusions Willis (1978,1980a)has considered the construction of bounds of HashinShtrikman type for a matrix containing inclusions. As in Section IV,B, the inclusions occupy regions V,, A = 1, 2, . . . , all of identical size and shape, and have moduli L1 while the matrix has moduli L2. It is convenient, however, to let Lo remain arbitrary and to let z* take the constant values zl and z2 in the inclusions and matrix, respectively. The term (z*, rz*) in F(z*) may be simplified by noting that the mean value of Tz* over V is zero, so that

this gives integrals only over the inclusions. With the infinite-body form (3.33) for r, it follows that

(4.77)

37

Overall Properties of Composites

+ c 2 t 2 .It is useful to define

where Z* = clzl

J

V,’

dx

J

dx’r(x - x’) = P

( A = B),

(4.78) ( A B), (so that P is given by (4.34) when the inclusions are ellipsoids). Then, Eq. (4.77) gives VA

+

VE

= QAB

(r*,rZ*)= (b/v)1( T i - 7 2 ) p ( z 1 - z*) + A

+ (71 - t2)J“-“A

dxf,(x’)(z2

1

(ti

BfA

-~ Z ) Q A B (~ ~7 2 ) (4.79)

- T*).

Replacing the right side of (4.79)by its expectation value and evaluating the other terms in F(t*) now gives F(z*) =

c cr[zr(Lr Lo)-%, 2

-

r= 1

x

{

czP+

sV

-2

4 + c1(z1 - z2)

~XBCPBIAQAB - clrA(xB)]]

(zl

-

t2)3

(4.80)

since t 1 - t* = c2(t1- t 2 )and t 2- T* = cl(tz - ti).z1 and z2 are now chosen to extremize (4.80). With the notation p’ = p

+

s

v-v,

the associated variational estimate

~XB[PB(AQAB - clfA(xB)],

(4.81)

for L takes the form

2

It can be verified that this is identical to (4.66), with PI = 2 and P replaced by P . For a composite containing more than one type of inclusion, an expression like (4.81)for F(t*)could be derived, but this, generally, would contain more than one tensor P and its extreme value would not yield such a simple expression for Z. Suppose, now, that P B I A is isotropic or, more generally, has “ellipsoidal” symmetry. The reasoning that gave Eq. (4.48) shows that the integral in (4.81) (for which the region VA is excluded) is zero, so that P = P. If, in addition, Lo is chosen as L,, elementary manipulation shows that (4.82) is identical to (4.54),which was obtained by use of the quasicrystalline approximation. Thus, use of this approximation, whose validity can be guaranteed only at low concentrations, in fact yields a variational estimate for L at any concentration, as does (4.59) for the n-phase material.

J . R. Willis

38

It is perhaps worth mentioning that the estimate (4.82),which was given by Willis (1980a),contains allowance for arbitrary inclusion shapes and distributions, through the tensor P'. In the particular case of "ellipsoidal" symmetry, when P = P, (4.82) [or, equivalently, (4.66)] yields bounds for the overall moduli of composites containing aligned ellipsoids. Examples for the related problem of conductivity (reviewed briefly in Section VI,A) were presented by Willis (1977), while for elasticity, estimates for L based upon (4.82)have been obtained by Willis (1980~)in the course of a study of wave propagation. Interactions between inclusions can also be treated from the variational principle. If, as in Section IV,B, the comparison material is identified with the matrix and the restriction of z* to VAis called zA,let (4.83) where the constant zo and the function fare to be determined. Substituting (4.83) into (4.67) gives a long expression with terms involving up to four inclusions, whose joint probability densities are unknown. However, an approximation correct to order cf is obtained by keeping terms involving only one or two inclusions. The resulting estimate for F(z*) has been given by Willis (1978); its extreme value, again correct to order c:, yields the estimate

z = L, + C,[(L,

- Lz)-l

+ PI-' + Cl[(L1 - Lz)-' + PI-'

x P[(L, - LJ-1 + P I - ' x Jv-,

+ C'[(L' X

Jv-",

- C,[(L,

- LJ'

+PI-'

~XBCPB~AQAB - CITA(XB)I[(LI- LJ-' +PI-' - L,)-'

+PI-'

~ J ' B IA Q A B [(L - LA-' I

x [(L, - La)-'

+ PI-'.

+ P + QABI-'QAB (4.84)

For a genera1 inclusion shape, (4.84) is not an exact estimate of L, even to order cl, because (4.83)takes z* constant over any inclusion. For ellipsoids, however, (4.84)agrees exactly with (4.45)to order c1 and so provides a variational estimate of the coefficient of cf in the expansion for L. Furthermore, this estimate of c: is a bound whenever (L, - L,) happens to be definite. Willis and Acton (1976) produced an approximation I,for a matrix containing an isotropic dispersion of spheres, by solving (4.47)by iteration. This amounts, in (4.84),to neglecting QAB in the term [(L, - L2)- + P + QAB]in the last integral. The variational estimate (4.84) was plotted out for an

Overall Properties of Composites

39

isotropic distribution of spheres by Willis (1980a). It compared quite well with the exact coefficient of c f computed for particular examples by Chen and Acrivos (1978b) and generally followed the estimate of Willis and Acton (1976), though the latter actually gave an estimate for it just less than that obtained from (4.84)in one case where this was, in fact, a strict lower bound. 3. Examples The various estimates and bounds so far given, although explicit, still usually require a small computation for application to a particular composite. The reader is referred to the original papers cited for such results. Here, a few limiting cases which may be discussed algebraically are presented for illustration. If the comparison material is taken isotropic, so that Lo = (3x0,2p0), and PSIror PBIAis isotropic as well, the Hashin-Shtrikman estimate (4.66) for L involves the tensor P, which is defined by (4.39),with K , p now taking the values K ~po, . For a composite comprising n isotropic phases with moduli L, = (3~,,2p,),(4.66) yields directly ?, = (3i?,2,ii),where (4.85) in which

The estimates R, ,iiprovide lower bounds for R, iiwhenever K~ < K, and p o < p, for all r, and upper bounds whenever K~ > K, and p o > p, for all r. The best lower bounds, correspondingly, are obtained by taking ico = K [ = min(ic,) and po = p l = min(p,); this is a limiting case in which the terms for which K , = K~ and p, = p l in (4.86) are replaced by zero. Likewise, the best upper bounds are those for which K~ = rcg = max(K,) and p o = pLB= max(p,). These bounds were given by Hashin and Strikman (1963) with the implicit restriction that ~ , , phad [ both to be obtained from the same phase, and K ~ pLssimilarly. The restriction was removed by Walpole (1966a), who also gave the bounds in alternative but equivalent forms. Another example for which the isotropic form of P is relevant is provided by an isotropic polycrystal. This is regarded as a limiting case of a multiphase material in which the moduli L, of the rth phase are obtained by a rotation, Q, say, of the moduli L, of a single crystal, taken relative to some convenient reference frame. The parameter r is now a continuous variable

,

J. R . Willis

40

which defines crystal orientation. All orientations have equal probability, which ensures overall isotropy of the polycrystal; polycrystals with texture would either have this restriction relaxed, or the restriction to isotropic PSI, relaxed, or both, with corresponding complications for the analysis. The tensor T, defined by (4.64)is anisotropic but its mean value (T) is isotropic, so that

(T)

= (37%

2pd,

(4.87)

say: with this definition, 2, ji are still given by (4.85).Heavy calculation is avoided by noting that I C ~p ,T can be obtained from the scalars ( T j j k k ) , (Tijij), as in (4.38). The corresponding scalars formed from T, are independent of crystal orientation. They are, therefore, identical to the required mean values and can be calculated for L, = L,. Consider, for illustration, a cubic crystal for which, relative to the cube axes, the constitutive relation takes the form Okk

= 3ucekk9

g11

- 622

= 2pc(ell - e22),

GI2

= 2p:e12,

(4.88)

with equations for the remaining components being obtained by cyclic permutation of the suffixes. It follows immediately that (Lc)iikk

= 9Kc5

(Lc)ijij = 3 K c

+ 4& +

(4.89)

The notation L, = ( ~ I c 2pc,2p3 ,, introduced by Walpole (1966b) allows products and inverses to be worked out directly so that, if T, represents T, referred to the cube axes, Eq. (4.64)gives

+ 3Kp}-l,

T, = ({[3(4 -

{ "4 - P0)I - + 211,)

{[2(pc - p0)l-l

+ 2pP }-' , (4.90)

-

Then, from (4.38)and (4.89), 3% = ([3(ICc -

+ 3Jcp}-1,

2PT = 5{[2(Pc - po11-l

+ 2pp}-l + ${[2(p:

- p(J1-l

+ 2pp}-1

(4.91)

and the estimates (4.85)are completely specified. It is easy to show that the tensor L = ( 3 ~2p, , 2p') is positive definite if and only if all of ic, p, and p' are positive. Bounds for R, p, therefore, follow by substituting into (4.91)the appropriate extreme values for ice, p o . The unique choice for x0 is IC,. This gives K~ = 0 for any p o and correspondingly, from the first of Eqs. (4.85), IC = K , exactly. The two possible extreme choices for p o are p, and pi: the greater provides an upper bound for p and the lesser a lower bound. When p o = pC,the second of Eqs. (4.85)gives, after rearrangement, p = pl, where

Pl = P c + 3{[5/(& - Pc)l + 8PP)

(4.92)

Overall Properties of Composites

41

while, when po = p:, it gives Ji = Ji, , where

p 2 = p: + 2{[51(~, - p:)i + 12~;)- l.

(4.93)

In (4.92),p p is evaluated with Lo = (3ic,, 2pc),while in (4.93),pb represents pp evaluated with Lo = ( 3 q , 2p3. The bounds (4.92), (4.93) were first given by Hashin and Shtrikman (1962b) and were discussed further by Walpole (1966b). Examples for which P is anisotropic are usually rather complicated, but two simple limiting cases are worthy of mention. Both involve spheroids with axes of symmetry parallel to the 3-axis, distributed with corresponding “spheroidal” symmetry. The tensor P is then transversely isotropic and it is helpful to introduce the notation of Walpole (1969), in which a transversely isotropic tensor A is written

b, b, d, 2f 29)

A=

if the relation a

= Ae

613 =

A-

(4.94)

between symmetric tensors G, e implies

3(6,, + 6 2 2 ) =

With this notation,

3

&ll

+ e22) + be337

2ge13,

‘ = (d/2A, - b/2A, - b/2A, a/A, 1/2f, 1/29)

(4.96)

where A = ad - b2. For a composite containing aligned “need1es”of length 21 and radius ~ 2 , with E << 1, in an isotropic matrix with Lame moduli 1, p, Willis (1980b) showed that

which may be substituted into (4.66) to yield estimates for E. The terms of order E’ are insignificant unless the terms of order unity cancel. This occurs when the needle-shaped inclusions are rigid and Lo is taken equal to L2,

J. R. Willis

42

the tensor of moduli for the matrix. This yields a lower bound for E (the upper bound being infinite) which is obtained from (4.54) with L, infinite:

r, = Lz + clP-l/(l

(4.98)

- c,).

If the needles are distributed at number density n,, their volume concentration is c1 = 4nn,Z3~2/3,which tends to zero if n, is taken finite. The tensor P-', however, contains terms of order E - and ~ gives

where 1, ,u here denote the Lame moduli that define L2. This is exactly the low-concentration estimate (4.45), but it provides a bound at any concentration. Similar degeneracy occurs for a matrix containing aligned platelets of a and thickness ~ U E where , E << 1. In this case, Willis (1980a) has shown that P

N

+

(0, 0, 0, (1 2 p , 0,1/2p) x

+ P V + 2P) IIE

($(A + 3 / 4 9 -$(A + PI, 4 1 + PI, $(A - P), &(31 + 7P), Q(31+ 4 d ) .

(4.100)

If the platelets are rigid, a lower bound estimate for the tensor of moduli of the composite is

-

2 L2 + (32/3)n,a3p(ll + 2p)([3(1 + 3p)]-l,

O,O, 0,2/(31

+ 7p),0).

(4.101)

For penny-shaped cracks, regarded as vanishingly thin spheroidal cavities, Eq. (4.54) produces a lower bound estimate for L which takes a particularly simple form if it is inverted to give an upper bound M for the tensor of compliances M. This is

M - M 2 + 4(n1a3)(a ") (0,0, 0, (i + p)-', 0,2/(3R + 4p)), (4.102) 3P which is linear in n1u3,although f, is not. The results (4.99), (4.101), and (4.102) were given by Willis (1980~). +

V. Self-Consistent Estimates Consider a composite comprising a matrix with moduli L,,, in which are embedded n different types of inclusions.Phase r is distributed at volume concentration c, and has moduli L,. The inclusions that comprise phase r

Overall Properties of Composites

43

have identical shape, except for r = n + 1 which defines the matrix. Suppose that the composite is subjected to a uniform mean strain Z which generates the mean stress B = LZ.If the average values of the stress and strain over the rth phase are defined as b,, e,, it follows that br

(5.1)

= Lrer

r=l n+ 1

ii =

1 c,L,e,.

r= 1

(5.3)

Eliminating cn+,en+ gives n

s=Ln+,s+

1 cr(Lr-Ln+l)er

r=l

(5.4)

so that, if e, = A,V,

(5.5)

the tensor of overall moduli is given by n

z=Ln+1+

1 cr(Lr-Ln+1)Ar,

r= 1

(5.6)

exactly. Sections IV,B and IV,C give, in effect, estimates for the tensors A,. In particular, z, might be calculated by solving Eqs. (4.75) and then A, deduced from the relation e, = (L, - Lo)-%,. A simpler approximate procedure is to estimate e, as the mean strain in an inclusion of rth type, embedded in an infinite homogeneous matrix with moduli equal to the overall moduli L of the composite. If, in particular, the rth inclusion is an ellipsoid, Eq. (4.31), with L1, L2 replaced by L,, L, respectively, gives z, = [(L, - L)-1

and, correspondingly, A, = [I

+ Pr]-lZ

(5.7)

+ P,(L, - L)]-’,

where P, represents the tensor P, defined by (4.34),for the rth ellipsoid. For inclusion shapes other than ellipsoids, the solution of (4.24) is required. Formally, it is convenient to define z, to be the mean of z over the inclusion and to define P, so that (5.7) is true. With this agreement, substituting (5.8) into (5.6) yields an equation which must be solved for L because it appears also on the right side. This “self-consistent” prescription for estimating L is essentially that given by Hill (1965b) and Budiansky (1965).

J. R. Willis

44

Suppose, however, that the composite has no clearly defined matrix phase to be singled out for special treatment or, equivalently, that c , + ~= 0. If the shapes of the phases can still be distinguished, it remains possible to estimate r, and A, by (5.7) and (5.8), but (5.6) no longer applies. Instead, a variety of desirable restrictions may be identified. First, Eq. (5.2) (with c,+ = 0) requires that n

1 cr[I + Pr(Lr- E)] -

r= 1

= I.

(5.9)

Also, a = LO coupled with (5.3) implies

c c,L,[I + P,(L, - E ) - y n

r= 1

and finally, C = Lo%!

(5.10)

=E

+ ( r ) , with Lo = E, implies c c,[(Lr - E)-1 + PI]-' = 0. n

(5.11)

rz1

Equations (5.9)-(5.11) provide three alternative methods for estimating L. They are not the same, in general, for Eq. (5.9) may be expressed in the alternative form n

1 crPr[(Lr-

r= 1

L)-1

+ Pr]-' = 0

(5.12)

and Eq. (5.11) is equivalent to

cs[I

+ P,(L, - L)]y

I-

= L.

(5.13)

If, however, Pr = P for all r, Eqs. (5.9)-(5.11) are equivalent. Furthermore, in this case, if cn+ # 0, (5.6) can also be placed in any of the forms (5.9)-(5.11) with n replaced by n + 1, even though the mean strain in the matrix is not estimated from the solution of an inclusion problem. A different perspective on the self-consistent method has been given by Willis (1977). Instead of estimating the tensors P, by solving inclusions problems, the Hashin-Shtrikman variational principle provides the optimal piecewise-constant polarization z* for estimating L, as the solutions r,, r = 1,2,. . . , n (or n + 1 if there is also a matrix phase) of (4.75), for any chosen comparison material with moduli Lo. If these equations give zr = s,e,

(5.14)

Overall Properties of Composites

45

the associated estimate for L is ULO) = Lo + (S),

(5.15)

where (S) = c c , S , . An estimate of self-consistent type for f, now follows from the assumption that f,(Lo) will be closest to E when the comparison material is identical with the “overall” material, that is,

Z(L) = E,

(5.16)

(S) = 0.

(5.17)

or, equivalently, Equation (5.17)states that the mean polarization is zero. This is also implied by (5.11) although, in that equation, the polarization is estimated differently, so that (5.11) and (5.17) are not equivalent. It may be noted that Eq. (5.17) contains information on the two-point probabilities P,, through the solution of (4.75) and requires no commitment to particular inclusion shapes. Suppose now that all of the probabilities P,, have the same “ellipsoidal” symmetry. The estimate (5.15) then reduces to (4.66) and (5.17) implies that (T) = 0. This is (5.11) precisely, with P, = P for all r. Thus, for such composites, the variational approach leads to any one of the equivalent prescriptions (5.9)-(5.11) without requiring any assumption for the shapes of inclusions. Naturally a problem remains if different inclusion shapes are known to occur. There is then little option but to accept (5.6) with A, estimated by (5.8) although, in principle, a generalization of the variational method given in Section IV,C,2 could be developed, if the relevant conditional probability densities were known. It may be remarked finally that, for the composite treated explicitly in Section IV,C,2, namely, a matrix with a single population of inclusions, the self-consistent equation (5.16), with t(Lo) given by (4.82), leads to any of (5.9)-(5.11) with P, replaced by P , which need not relate to “ellipsoidal” symmetry. Two simple examples of the use of the self-consistent method are now presented, for illustration. First, for a matrix containing spheres, distributed isotropically, L is estimated as the solution of (5.9), (5.10) or (5.11), with P given by (4.32),in general, or by (4.39)if L is isotropic as will now be assumed. If there are n - 1 different inclusion types, all isotropic with L, = (3rc,, 2pJ and the matrix ( r = n) is likewise isotropic, the prescription (5.11) implies that rcT and p T , defined by (4.86),are both zero when rc0, p o are chosen as the “overall” moduli. For a matrix containing a single family of inclusions, so that n = 2, elimination of il from the self-consistent equations leads to a quartic equation for p. This will not be given, but, in the particular case where L1 = 0, so that the inclusions are actually cavities, it can be shown

J . R. Willis

46 that

(5.18) while p satisfies the quadratic equation

16p2 + p[Kc,(3

- c1)

- 4p2(4 - %I)]

- 3pzK2(1 - 2 ~ 1 = ) 0.

(5.19)

In the limit xZ -, 00, corresponding to incompressibility of the matrix, these equations imply 4(1 - ~ 1 ) (-l 2 ~ 1 ) 3(1 - 2 ~ 1 ) (5.20) Pz, F = 3 - c1 PZ > c1(3- c1) as given by Budiansky (1965). Equations (5.20) predict that the overall moduli reduce to zero at c1 = and indicate the need for some caution in applying the self-consistent method at high concentrations of inclusions with such extreme properties. It may be noted, for comparison, that the Hashin-Shtrikman estimates (4.85) give upper bounds for IC, p when Lo is chosen identical with Lz . When the matrix is incompressible, these reduce to (5.21) = 4pz(1 - c1)/3c1, p = 3/lz(1 - C,)/(3 + 2Cl). R=

Lower bounds would require Lo = L1 = 0, giving E = 0. As a second example, a self-consistent estimate is obtained for the overall moduli of the polycrystal for which bounds were given in Section IV,C,3. If the associated single crystal is cubic, with L, = (3rcC,2pc,2p3, the selfconsistent prescription (5.17) requires that and pT are zero when K ~po, take the values R, p, where p T are given by (4.91).The equation for I C gives R = K , and elementary manipulation of the equation p T = 0 then gives the cubic equation

8p3

+ ( 9 ~+, 4p,)ji2 - 3pi.i~~ + 4pC)F- 6~,p,p;= 0

(5.22)

for p. This result was first given by Hershey (1954) who assumed spherical grains and applied the prescription (5.10).This led to a quartic equation for p given, in present notation, by (5.22)multiplied by the factor ( 8 p + 9q). It is perhaps worth emphasizing again that (5.22)follows from the variational approach directly from the assumption of isotropy of the polycrystal, without reference to grain shape. The self-consistent method has been applied to a variety of examples less straightforward than those described above. Kneer (1965) studied textured polycrystals, in which grains were modeled as spheres, but not all crystal orientations were equally likely, so that L was anisotropic and P was computed from (4.32).Walpole (1969) considered composites containing transversely isotropic “needles” and discs whose properties were such that end

~

Overall Properties of Composites

47

or edge effects were not significant and P could be estimated asymptotically from one- or two-dimensional problems. Budiansky and OConnell (1976) studied a body containing flat cracks with elliptical boundaries, oriented randomly so that E is isotropic. Circular cracks oriented to make E transverselyisotropic were consideredby Hoenig(1979),andLaws and McLaughlin (1979) studied composites containing short fibers, modeled as aligned spheroids; here again, P was computed for a transversely isotropic matrix. VI. Generalizations The discussion so far has been expressed in terms of elasticity theory. However, the methods (and, in many cases, even the results) of the preceding sections apply to a variety of mathematically similar, though physically distinct, problems. First, if the phases are taken as incompressible and u is interpreted as velocity rather than displacement, the original problem becomes that of finding the overall viscosities of a fluid suspension in the limit of Stokes flow. Of course, for the fluid problem, the configuration of the suspension evolves with the flow and the question studied here, namely that of finding the viscosity given the configuration, is by far the easiest aspect of the whole problem. Study of the evolution of the configuration is, not surprisingly, still at a primitive state of development; it is discussed in the review by Batchelor (1974), for example. The problem of determining the overall thermal conductivity of a composite is described by Eq. (2.16) and (2.17), if u is interpreted as the scalar temperature field, e is its vector gradient, n is the negative of the heat flux vector, and L is the second-order tensor of conductivities. In (2.17), which is now a scalar equation, f is the scalar heat sink field. The explicit form of the constitutive relation becomes, in place of (2.14), 0 I. =

L IJ. . eJ .= L..u IJ .J ..

(6.1)

The Green’s function G is a scalar and r is a second-order tensor operator whose kernel has components Tij(x,x’) = a2G(x,x’)/(dxiax>),

(6.2)

in place of (3.10). The infinite-body Green’s function G“ remains homogeneous of degree - 1 and every equation in which suffixes are not explicitly given applies directly to the conductivity problem. Willis (1977) has, for example, evaluated P for a spheroid in a transversely isotropic matrix and has calculated the conductivity of a body containing aligned spheroids. With different interpretations of n, e, u, and L, the thermal conductivity problem

48

J . R. Willis

is identical to those for electrical conductivities, dielectric constants and magnetic permeabilities. There are, in addition to these, a number of problems whose structure is similar to that of the elasticity problem but which require more explicit development; two such problems are discussed below.

A. VISCOELASTICITY The constitutive relation for a viscoelastic solid may be expressed in the form (2.16),so long as the tensor L is interpreted as the Stieltjes convolution operator

or, equivalently, in terms of generalized functions, L: e(t)+ J [ d ~ ( t t’)/dr]e(t’)dt’,

(6.4)

the function L(t) taking the value zero for t < 0. Either of (6.3) or (6.4) give, for the particular strain history e(t) = eoH(t), u = Le = L(t)eo,

(6.5)

showing that L(t) is the stress relaxation tensor. For a viscoelastic composite, therefore, the objective is to find the overall modulus operator L or, equivalently, the overall relaxation tensor L(t).Dually, if the operator L is invertible, with inverse M, the corresponding function M(t) is called the creep compliance tensor and the overall compliance operator M or creep compliance tensor M(t) could be sought. Formally, the viscoelastic problem can be reduced to the form of the elastic problem already considered, by application of the correspondence principle: for boundary conditions such as (2.19),with 5 now a function oft, application of either Fourier or Laplace transforms yields an “elastic” problem, in which the moduli are now functions of the transform variable. The elastic problem is solved and finally the transform is inverted. All of the methods given in Section III-V are applicable, therefore, so long as the relaxation tensor L(c) has the symmetry (6.6) which has been implicitly assumed for elastic bodies. The variational estimates of Sections III,B and IV,C are now merely estimates rather than bounds, however, because the transform of L is generally complex and cannot be taken as definite in the manner assumed for elasticity. Lijkl(t)

= Lklij(t)

Overall Properties of Composites

49

It is possible, in fact, to analyze the problem directly in the time domain. For fields 2, q defined over V , the appropriate generalization of (3.52)is the bilinear form

where * denotes the ordinary operation of convolution with respect to time. The definition (6.7) is time-dependent, but involves only values oft, q up to time t if they are both zero for times t < 0. Now define (6.8)

O(t)= t ( n ,e) = +(Le,e),

in analogy with (2.22).If e is associated with displacements that conform to the boundary condition (2.19), with i? now a function oft, it follows that (6.9)

o(t)= f ( i 3 , E ) = $(LE,E), as in (2.23). Also, the symmetry of the bilinear form

(Le,,e,) = 1/V

sv

[(dL/dt) * e,]

* e,dx

(6.10)

yields the variational principle for the displacement boundary value problem, that the functional (6.11)

8*(t)= +(Le*,e*),

where e* is the strain associated with any displacement taking the prescribed boundary values, is stationary at each t for the actual strain field, e. Variational principles of this type are summarized by Leitman and Fisher (1973). Following Laws and McLaughlin (1978), it is useful to define the Green’s function G(x, x‘, t) for a homogeneous comparison body with modulus operator Lo so that it satisfies

+ 6,6(x

- x’)H(t)= 0,

x E V,

t

> 0,

(6.12)

Gip(x,x’,t) = 0,

x E at‘,

t > 0,

(6.13)

Gip(x,x’,t) = 0,

x E V,

t I 0.

(6.14)

The representations (3.6) and (3.7) then hold for the solution u(x, t ) of the viscoelastic boundary value problem (3.2)and (3.3),so long as G is interpreted as an operator of the same type as L. With this preparation, the whole of Sections 111-V apply also to the viscoelastic case, except that variational

50

J. R . Willis

estimates no longer provide bounds. Of course, the difficult task is to evaluate the operators explicitly.In practice, this is likely to be performed using transforms. The most complete study to date is that of Laws and McLaughlin (1978), who evaluated the operator P for spherical and circular cylindrical elastic inclusions ( E glass) in a viscoelastic matrix (a cold setting epoxy polymer) and obtained corresponding self-consistent estimates for the overall creep compliance tensor.

B. THERMOELASTICITY Static problems of thermoelasticity concern, in addition to stress and strain fields a,e, heat flux and temperature fields q, 8 which are related as follows: u = Le - 18,

(6.15)

q = KV8.

(6.16)

In (6.15) I is the stress-temperature tensor and L is the tensor of isothermal elastic moduli. Temperature 8 is measured relative to a reference temperature TR,say, at which the body is stress-free when its strain is zero. The tensor K in (6.16)is the tensor of thermal conductivities and is positive definite. In the discussion at the beginning of Section VI, it was implicitly assumed that K was symmetric. This is the usual assumption, though Carlson (1972) has pointed out that it is not a requirement of thermodynamics. In the absence of heat sources, divq = 0 and this equation, together with (6.16) and the thermal boundary conditions, define a purely thermal problem which can be solved independently of conditions on stress or displacement. In particular, if 8 takes a prescribed constant value over aV, then it is constant over I/. Whether or not this condition is adopted, 8 may be regarded as known in (6.15) to generate, with (2.17) and the stress or displacement boundary conditions, a purely mechanical problem in which 8 provides one of the sources of stress. Restricting attention now to the case of uniform temperature 8, zero body forces and displacement boundary conditions, the principle of minimum energy states that the actual displacement field is that which minimizes the Helmholtz free energy of the body, whose density +I! is given by

i,b = +eLe - 8k - iff?’.

(6.17)

The scalar f is introduced in (6.17) to allow a complete thermodynamic description. Equation (6.15) is equivalent to a=

(6.18)

51

Overall Properties of Composites

and, if q denotes entropy density, q=

-a$/as

= Ie

+ fe.

(6.19)

The specific heat at constant strain s(e) is then s(e) = T,(aq/ae), = ~ Dually, with the definitions g =f

m = MI, Eqs. (6.15) and (6.19) may be rewritten

e = Ma q=

~

f

.

+ Im,

(6.20) (6.21)

+ em,

(6.22)

+ gf?,

(6.23)

and the specific heat at constant stress s(u) is s(u) = TR(dq/dO),

=

(6.24)

TRg.

The tensor m is the thermal expansion tensor. The only bounds so far published for thermoelastic properties of general composites have been obtained by substituting into the minimum energy principle a uniform strain field and, dually, substituting into the complementary energy principle a uniform stress field (Schapery, 1968; Rosen and Hashin, 1970), generalizing the procedure given in Section II,B,l where the Voigt and Reuss bounds for L were derived. This work will not be reviewed, but rather, an outline of generalizations of the methods of Sections III-V will be given. A uniform comparison material, with properties Lo, I,, fo is introduced, relative to which the constitutive relation (6.15) may be written u = L,e

- 1,s

+ Z,

(6.25)

where the stress polarization z satisfies

z = (L - L,)e - (1 - lo)&

(6.26)

generalizing (3.24). If 6 is uniform and the boundary condition (2.19) is imposed, the strain e may be represented, analogously to (3.8), in the form e = e - Tr,

(6.27)

to which perturbation theory may be applied, as in Section III,A, if L - Lo is small. It is more interesting, however, to eliminate e from (6.26) and (6.27) to give equation [(L - L J - ~+ rlt = E - (L - L,)-'(I

- lop,

(6.28)

J. R . Willis

52

which is (4.1) with a particular e0 on the right side. The Hashin-Shtrikman principle, therefore, shows that the functional ~ ( z * )= (z*, ( L - L,)-'Z*)

-2(r*, E )

+ (z*,rz*)

+ 2(z*, (L - LJl(1

- I,)@

(6.29)

is stationary when z* = z, the actual polarization and, furthermore, that F(z) is a minimum if L - Lo is positive definite or a maximum if L - Lo is negative definite at each point of V . The extreme value F(z) is expressible in terms of the mean $ of the Helmholtz free energy. In fact, ~ ( z= ) -(0,q

+ (2,(L - L,) - ](I - io)e),

(6.30)

which gives, upon substituting for z using (6.25) in the first term and (6.26) in the second,

F(z) = PLOP- 281,C - (a,E)+ (e,le) - ((I

- I,)(L - L , ) - ~ I- io))e2.

(6.31)

Now from (6.17) 2$ = (a,4 - (e,W - ( f > 0 2 ,

(6.32)

where (f ) denotes the mean off over V. The stress a is divergence-free and so (o,e) = (u,E)= (a,E).Therefore, ~ ( 7= )

2+,

+

foe2

- 2~

(6.33) - [(f> + ((1 - lo), (L - Lo1-V - 10))]e2. Estimates for overall properties L, i,7 follow if $ is defined to have the

form (6.17). With a temporary distortion of the notation used in Section IV,C,l, for an n-phase material, let z, denote the mean of z over the phase r and let a,, e, denote the corresponding means of stress and strain. Then, by averaging (6.25) over phase r, a, = LOer- 1,8 + 7,. (6.34) and, averaging (6.26) and rearranging, e, = (L, - L,)-'[z,

+ (1, - l,)O].

(6.35)

Hence, by substituting into (6.32), 2iJ = ZLOE - 281,s -

+ &r,s

ecc,(i, - i,ML, - L,)-

it,

- C c r [ f i + (1, - 1o)tLr - LJ-lOr

- i,)]e2.

(6.36)

Overall Properties of Composites

53

Now z, is a linear function of Z and 8, say (6.37)

zr = S,F - s,e.

If 3 has the form (6.17), it follows then that

-

+ CCrSr, + &[s, + (1, - IO)(L, - Lo)-'S,],

L = Lo

I = lo

9 = C C r C f i + (1,

- Lo)-'(lr - 10 - s r ) ] .

-

The average of (6.34) gives

a = (L, + ~ C , S , ) F - (I, It is desirable, therefore, if a = a$/%, that Ccrsr

= Ccr(Ir -

(6.38) (6.39) (6.40)

+ Cc,s,)e.

(6.41)

- LO)-'sr*

(6.42)

This relation is in fact true: it can be shown to follow from the selfadjointness of (6.28). The relations (6.38)-(6.42) have been given, in terms of stress concentration factors which relate e, to z! and 8, by Laws (1973), who deduced the symmetry (6.42) from the minimum energy principle. The formulas of Laws may be related to those above by expressing e, in the form

e, = A$

- a,&

(6.43)

Equations (6.35) and (6.37) imply

A, = (L, - L0)-'S,,

a, = (L, - Lo)-'(sr + lo - I,),

(6.44)

and (6.42) is equivalent to CcJrAr = C c r ( I r

+ La,).

(6.45)

This follows upon use of the relations Cc,A, = I,

'&a,

= 0,

(6.46)

which are implied by averaging (6.28) over V . If z, is now estimated by substituting z*(x) = C z r f i ( x )

(4.72)

into F(z*) and extremizing, as in Section IV,C,l, a result of the form (6.37) is obtained and the symmetry of the algebraic equations [which are just (4.75) with more complicated right side] implies that the approximate S,, s, so generated satisfy (6.42). Also, the extreme value for F(z*) is given by (6.30), except that z, are now those obtained via (4.72). Completing the

54

J . R . Willis

algebra gives

z(L - L)E - 200 - i)z - ( f - 7 ) I~ o

(20)

z,

(6.47)

whenever (L, - Lo) is positive (negative) definite for all r, where l, and are approximations that correspond exactly to (6.38)-(6.40) with S,, s, now defined as the Hashin-Shtrikman approximations. The inequality (6.47)yields the Hashin-Shtrikman bounds for L that were discussed in Section IV,C,l and also gives directly corresponding bounds for and hence for the specificheat at constant strain. Bounds for 1are obtained indirectly, by considering the quadratic form. To illustrate this, suppose that the composite (though not necessarily its constituent phases) is isotropic, so thatland interact with i?only through its dilatational component. Definiteness of the inequality (6.47) then requires that

f’

7

( B - ZSY I (R - ?s)(7 - TI,

where R, iz represent the bulk moduli associated with E, L and B, diagonal entries of ‘i;7. The inequality (6.48)may be written - [(R - E)(f -

Jo]l’z I 7J I B + [(R

- E)(f - f)]l’z

(6.48) are the (6.49)

and bounds not involving E,f (whose values are not known) may be obtained by replacing R - iz, f - f by the differences between the Hashin-Shtrikman bounds for i?, f,respectively. The bounds produced from the classical energy principles by Schapery (1968) and Rosen and Hashin (1970) also required considering a quadratic form but seemed always to involve the exact overall moduli L. It should perhaps be mentioned that the full development of (6.47) is unnecessary for a composite with only two phases. Laws (1973) derived the exact result

-

1 = (L - LJ(L1

- LZ)-’I1

+ (L - L,)(LZ - Ll)-’lz,

(6.50)

by eliminating the tensors A,, a, (or, equivalently, S,, s,) between equations (6.38),(6.39),(6.42), and (6.46), so that a bound for E immediately induces a bound for 1.It may be noted that, since the Hashin-Shtrikman estimates for S,, s, conform exactly to (6.45), (6.46), the estimate corresponding to (6.39) for is expressible in the form of Eq. (6.50), with L replaced by E, so that this estimate provides bounds automatically. Finally, Laws (1973) proposed the estimation of A,, a, (or, equivalently, S,, s,) by the self-consistent method of embedding an inclusion of rth type in a matrix whose properties are those of the composite, as described for the mechanical problem in Section V. The remarks of that section are applicable here : the usual self-consistent method has a variational interpretation for composites with “ellipsoidal” symmetry and, for more general

Overall Properties of Composites

55

composites, a “variational” prescription, which allows explicitly for the statistics through the probability P,, ,is available as an alternative. Explicit self-consistent formulas for a composite containing spherical inclusions were given by Budiansky (1970). VII. Problems Which Lack Convergence Beginning in Section III,A, one of the major points of principle to receive emphasis has been the desirability of formulating a boundary value problem exactly for a finite body and only afterward making any simplification that is allowed when the body is, in fact, large in comparison with microstructural dimensions. This led, for example, to the prescription (3.33) for the limiting form of the operator r which is not just a direct replacement of r by the operator I‘” derived from the infinite-body Green’s function. The kernel P ( x ) decays like 1x1- at large 1x1 and, if applied directly to z, would lead to an integral that is only conditionally convergent. This type of complication is altogether more severe for problems which generate operators whose infinite-body forms decay more slowly at large 1x1. It is the writer’s view that precise formulation of such boundary value problems prior to the development of infinite-body approximations is highly desirable, at the very least. Such initial formulations have, however, been very much in the background if they have been acknowledged at all in most studies carried out to date, both for problems of “overall modulus” type and for the problems that will now be described. The problem with a “badly behaved” Green’s function that has received most attention is that of determining the resistance to flow of viscous fluid through a fixed bed of obstacles. Macroscopically, the resistance takes the form of a body force f related linearly to mean fluid velocity (u). If the fixed bed is statistically isotropic, the relation reduces to

f = -paZ(u),

(7.1)

where p represents the viscosity of the fluid and the Darcy coefficient is to be estimated. A fixed bed of spheres has been considered by several authors, starting with Brinkman (1947) who performed a self-consistent calculation, rather like that described in Section V. Brinkman’s idea has been developed further by Lundgren (1972) and Howells (1974) to make explicit allowance for correlations between sphere positions. Low-concentration approximations for the Darcy coefficient are of considerable theoretical interest (largely because of the convergence difkdty), though probably of limited practical

56

J. R. Willis

value; they have been considered by Childress (1972), Hinch (1977), and others. The convergence difficulty may be explained by discussing the flow of viscous fluid past a single bounded obstacle. Far from the obstacle, the induced perturbation to the mean flow has the “Stokeslet” pattern generated by an isolated body force. At distance r from the obstacle, this is of order r - ’ . Therefore, if a second obstacle is placed in the flow, its interaction with the first is of order r - l if their separation is r, and any attempt to allow for a distribution of obstacles by summing interactions between pairs is doomed to failure through divergence of the integral. The problem is even worse in two dimensions: Green’s function is of order lnr and the problem of an isolated cylinder perturbing a uniform mean flow cannot even be defined. Self-consistent calculations and perturbation expansions are the only type that have been applied so far to the fixed-bed problem. It is likely, however, that a variational approach could be developed: all three methods have been applied by Talbot and Willis (1980) to a problem with similar mathematical structure, which will now be outlined. A. A MODELPROBLEM:DIFFUSION TO A RANDOM ARRAYOF VOIDS Irradiation of metal by neutrons, electrons, or ions generally produces overall distortion. This is often associated with the atoms rearranging themselves so that the specimen contains voids or gas bubbles, with consequent volume change. The full range of phenomena is complex and cannot be properly described here; the review of Bullough and Hayns (1978) provides background, together with details of a variety of calculations of selfconsistent type. It is apparent, however, that the basic rearrangement takes place by diffusion in the presence of various types of sink distribution, including the distribution of the evolving set of voids or gas bubbles. Talbot and Willis (1980) studied the simple model problem of determining the mean sink strength of a random array of voids in the presence of a single population of diffusing defects (gas atoms, say), when no other sinks were present. The material was taken to occupy the region V and was regarded as a matrix containing identical spherical voids occupying V, with centers xA, A = 1, 2,. . . , as described in Section II,A. Overall, for such a medium, if the defects have concentration c(x, t ) and are introduced throughout V at rate &x, t), the concentration is expected to be described by the diffusion equation aclat = b(v2c- k2c) + I?,

xE

v,

(7.2)

in which b is an overall diffusion coefficient and the sink term -bk2c, which is analogous to the Darcy resistance (7. l), represents the mean effect

Overall Properties of Composites

57

of the voids. If (7.2) is indeed a fair description, then it must apply in particular to the steady-state situation in which & is independent of x and t and no flux is admitted across dV, so that &/an = 0 there. The solution of (7.2) in this case is c = T, a constant, related to I? and k2 by Dk2C = I?.

(7.3)

The simplest way to estimate bk2is by the self-consistent method of considering the flux into a single void of radius a embedded in the “overall” material. If the boundary condition c = 0 is adopted at the void surface, this leads to the solution c = i ~ {1

- (a/r)exp[ - k(r - a)]},

(7.4)

where T is given by (7.3) and r represents distance from the void center. It may be noted at this stage that, if k = 0, (7.4) gives a perturbation of order r - while introduction of the overall sink term gives exponential decay and so would “screen” distant sinks from one another. The flux into the void is now calculated as

F = 4 n ~ (+ l k~)bF,

(7.5)

while (7.2) gives the sink strength per unit volume as bk2T. If the voids are distributed at number density n l , therefore, self-consistency requires that

k 2 = 4nanl(l

+ ka).

(7.6)

The overall diffusion coefficient b is not estimated so easily and it is usually identified with the diffusion coefficient D of the matrix. In contrast with the above, the steady-state problem is described exactly by the equations

DV2c + K = 0, dc/dn = 0, CEO,

x E V’ x E av, X € a v ~ ,

(7.7) A = 1 , 2,...,

where V’ denotes the part of V occupied by the matrix and the flux term K is related to R by

Iz = K(1 - CJ,

(7.8)

where c1 = 4nn,a3/3 represents the volume concentration of voids. The object is to solve the system (7.7) to obtain E, the mean of c over V (with c = 0 over VA),in terms of which (7.3) provides a precise definition of the sink term Bk2, analogous to the definition (2.21) of overall moduli.

58

J. R . Willis B. AN INTEGRAL EQUATION AND

PERTURBATION

THEORY

The problem defined by Eq. (7.7) may be expressed in the form of an integral equation, analogous to (4.1), if a Green’s function for the region V is defined so that V 2 G ( x , x ‘ ) = 6 ( x - x’) - l/V,

x E V,

aG/an= 0,

x E av.

(7.9) (7.10)

The presence of the distributed source should be noted in (7.9): although it is small when V is large, it is essential for the consistency of (7.9), (7.10). The solution of these equations is made unique by requiring that JV

(7.1 1)

G ( x , x’) dx = 0,

which also ensures that G ( x ,x’) = G(x’, x). Application of Gauss’ theorem in the usual way yields c(x’) = F

+ SaVv G ( x ,x’)q(x)ds - K’

JV,

G ( x , x’) dx,

(7.12)

where V, denotes the union of the regions V,, (7.13)

K‘ = K / D and q ( x ) = ac/an,

x Ea

v,.

(7.14)

The boundary condition c = 0 at the void surfaces now generates an integral equation for 4 ( x ) in which the additional unknown z is counterbalanced by the consistency condition

J,

q(x)ds = V’K’ = V(l - c,)K‘.

(7.15)

More conveniently, the integral equation may be viewed as an equation for q ( x ) given F, in which K is defined by (7.15). Proceeding as in Section IV,B, the restriction of q(x)to aVAis called qA(x).The integral equation may then be written in the form CGqA+F=O, A

X’EaVA,

A = l , 2, . . . ,

(7.16)

where operator G is defined by G : qA -+

LVA

+

dsG(x,x’)qA(x) K‘

VA

dxG(x,x’).

(7.17)

Use was made, in deriving (7.16), of (7.11) to replace the integral over V by integrals over V,, A = 1,2, . . . . Taking expectations of (7.16), keeping

59

Overall Properties of Composites

’V fixed, now gives

+

~ 4 2J ~ ~ B P B I , G ~ : B

+ z = 0,

XI

E

av,.

The term K‘ in the definition (7.17) of G now becomes
=[

~ ( l -CAI-’

{L~~+ qj(x)ds

L~~

J~~BPBIA

but, in the limit of large V , this is equivalent to

(K‘)= [v(l - ~ 1 > ] - ’JdxBPB

(7.18)

I

~ B A B w ~ X (7.19)

LvB&(x)ds.

(7.20)

The simplest approach to (7.18) is to apply the quasicrystalline approximation (4.51), which reduces it to an equation for q;. If it is assumed in addition that the operator G defined by (7.17) is translation invariant in the limit of large V , it follows that 4: is insensitive to translations of x A . Finally, if PBIAis isotropic, 42 reduces to a constant, qo say, and (7.20) becomes 3c1q0 ( K ‘ ) = 4na2n1q01 - c1 a(1 - cl)‘

(7.21)

The problem is now to find the appropriate representation for G. When q2 is constant over aV,, Talbot and Willis (1980) simplified the integrals in (7.17) by application of the mean value theorem for harmonic functions. The result is expressible in the form

1 - c1

-a(l

C1 + ic,) + Ix‘ - xA12+ 4na2G,(x,,x‘) 2a

(7.22) where GA(xA,x’)is a function that is essentially an “image” term and is small when V is large and x’ and xA are not close to dV. Substituting into (7.18) now gives, when I/ is large, 1 -c1 +Z=O,

X E ~ V ’ , (7.23)

J. R. Willis

60

in which the term involving V-’ cannot be ignored because it is integrated over I/. The left side of (7.23) is expected to be independent of x’: accepting that this is so, (7.23) may be averaged over dV,, which has the effect of replacing G(x,, x’) by G(x,, xA)- a2/6V. Finally, a convergent integral involving G is obtained by use of (7.11); the result is

{

[qo/(l - cl)] u - 4na2

s

dxB(PBIA

- ~,)G(x,,x,)

I

+ U C , (-~ ~ , / 5 ) = Z.

(7.24) If the voids have no long-range order, only points X~ close to xAcontribute significantly to the remaining integral and so G may be replaced by the infinite-body Green’s function. Thus, finally,

The author must confess to finding the above argument not entirely convincing. While it is true that the quasicrystalline approximation is likely to produce a result correct only to lowest order, the contribution of order c1 from the integral of a term of order V-l hardly lends confidence to the assumption of translation invariance on which the derivation depends. The estimate (7.25) nevertheless does have an honest interpretation, as will be shown below. Before proceeding, however, it may be remarked that interactions can be allowed for, as in Section IV,B,2, by taking expectations of (7.16) with two spheres fixed and making the closure assumption (4.55). A calculation equivalent to this was performed by Brailsford (1976) (who did not consider terms of order V-’, however). The “fixed-bed’’ problem has been approached in a similar way by Hinch (1977). In the present problem, an integral equation for (q:B - q z ) is generated, whose solution shows the same type of exponential decay as (7.4) and yields an estimate for 8 k 2 that agrees with (7.6) to order (c1)’l2. C. VARIATIONAL FORMULATION

Consider, in place of (7.16), the integral equation

syI

dxf(x)G(x,x’)

+ g(x’)= 0,

X‘ E

Vl,

(7.26)

where G is the Green’s function for V and V, is a subset of V . Equation (7.26) is self-adjoint, from the symmetry of G, and it is easy to prove that (7.27)

Overall Properties of Composites

61

for any function f* defined over V,. Equation (7.26) therefore implies the maximum principle,

(7.28) with equality when f* = f . Now Eq. (7.12) defines c(x‘) as a harmonic function within V,, which is zero because it is zero on the boundary aK. Equation (7.16) therefore applies for x E VA,A = 1,2, . . . and can be considered as a limiting case of (7.26). The corresponding limiting form of (7.28) is

{LVAds’qA(x’)+ K*

VA

dx’

B+A

(7.29)

where qA here denotes the restriction to aVAof any approximation q* to q and K* is the corresponding estimate of K’, as in (7.15). Any approximation q* thus bounds K , and hence bk2,from below. If the functions qA are now taken constant over aVA, the functional in (7.29)can be simplified using (7.22).The result is

a’ 15

- - (5 - c ~ ) ( K *+) ~2K*Z 5 K‘Z,

(7.30)

where (7.31) G A denotes G,(xA,xA) and G A B denotes G(xA,xB). The inequality (7.30) applies for any constants qA: they need not be all the same, as was assumed in the preceding section. It is essential, in fact, that they should vary with A for a useful result to be obtained.

62

J. R. Willis

Although the distribution of the voids is assumed statistically uniform (except close to dV), the Green's function is sufficiently unpleasant to discourage the replacement of the spatial mean over A in (7.30) by an ensemble mean over some chosen V,. Still, assuming that, for large V , the left side of (7.30) is independent of the sample CI,it is now replaced by its expectation value. If qAis taken to depend only upon the location x, of V,, the expectation involves integrals of probability densities for the location of up to four inclusions, because K* is defined by a sum. The object is to maximize the expectation of the left side of (7.30)by setting its variation with respect to qR equal to zero. The same result is obtained, in fact, by differentiating the left side of (7.30)with respect to qR and taking the expectation value of the result, conditional upon x R being fixed. With the approximation, valid when V is large, that the conditional expectation of any mean value (that is, V - times a summation) is equal to the unconditional expectation, Talbot and Willis (1980) showed that

''

ha3 -k 3V(1 - cl)

where

a2(5 - cl)

JdXApA4'- 15(1 -

cl)

-

(K*)

+- 0, 1 - c1 C

(7.32)

(7.33) Equation (7.32) shows that 4Ris independent of x R . This allows evaluation of the integral, to simplify (7.32) to

4R- a2(5 - Cl)(K*)/lS

+ I = 0.

Then, with the definition qA = [(l - c,)/3cl]a(K*)

+ 8,

(7.34) (7.35)

where (7.36)

J d X , P A 8 = 0,

Talbot and Willis (1980) showed that (7.34) requires

-

(a/3cl)(K*)(PB - nl) + eBPB 0

(7.37)

except, perhaps, in a layer close to aV where P, may deviate significantly from n, ; in this layer, the left side of (7.37) will have zero mean. Elsewhere, PB n, and so eB 0 and qA is constant, as assumed in Section IV,B. N

-

Overall Properties of Composites

63

Further analysis is required, however, because these statements are true at any chosen point x A , but terms of order V - may still contribute significantly to the functional. The end result is that

which is precisely consistent with (7.25).Since the expectation of the left side of (7.29) is quadratic, its extreme value is just one half of the linear term, which is ( K * ) l precisely. Hence, from the definition (7.3) of Bk',

bk' 2 (1 - c,)D(K*)R.

(7.39)

Assuming now that PslRis isotropic, it is convenient to express it in the form PBlR

where x = I

x -~ x,1/2a.

[

= nidx),

(7.40)

Then, g(x) = 0 for x 5 1 and (7.38) may be written

a2(K*)/3Z = 1 - SC, - c:/5

+ 12cl

J:

I'

(g(x) - 1)xdx

. (7.41)

The simplest possible choice for g(x)is g(x) = 1:this defines the "well-stirred" approximation, used in the low-concentration limit by Batchelor and Green (1972),Willis and Acton (1976), and others. Unfortunately, the approximation fails at high concentrations, predicting that (K*) becomes negative after passing through a singularity near c1 = 0.2. This is inconceivable, since it can only have resulted from the expectation value of the negative definite quadratic functional (7.27) having become positive. Equation (7.41) thus indicates that g(x) must necessarily rise above 1 when c1 is large. This is plausible at high concentrations, since knowledge that a sphere is located at some point must almost guarantee the presence of spheres centered 2a away. It is also borne out by more realistic models discussed by Talbot and Willis (1980). The most interesting is that derived from the Perkus-Yevick approximation to the pair function g(x)for a statistical-mechanical ensemble of hard spheres. The integral equation for g(x) produced by Perkus and Yevick (1958) was solved by Laplace transforms by Wertheim (1963). The transform requires numerical inversion (Throop and Bearman, 1965), but the integral in (7.41)can be deduced analytically from the small-argument expansion of the transform. With the Perkus-Yevick form for g(x), (7.41) remains finite, but, interestingly, the self-consistent estimate (7.6) lies below the lower bound when c1 > 0.2, if b is identified with D. This contrasts with the results of Sections IV and V which show that the self-consistent estimate of L always lies between the Hashin-Shtrikman bounds.

64

J . R. Willis

VIII. Wave Propagation Considerable effort has been applied to the study of wave propagation in random media: see, for example, the books by Chernov (1960), Uscinski (1977), and Ishimaru (1978). The vast majority of studies have been, however, for electromagnetic or acoustic waves, both of which are described essentially by the scalar wave equation. Methods developed for these fields of application can be carried over to elastodynamics at the expense of algebraic complication, but they need not always be relevant. Whereas, in the electromagnetic case, interest frequently centers on waves whose wavelength is short in comparison with distances over which variations in the properties of the medium are significant, the opposite limit is often the one that is required for elastic wave propagation through a composite. For example, a wave traveling with a speed 3 km s-’ typical of metals or rocks has a wavelength of 3 mm when its frequency is 1 MHz. Also, strong discontinuities in properties across phase boundaries are the rule rather than the exception in solids and limit the utility of methods that assume smooth variations. In any stochastic problem involving wave propagation, the approach invariably adopted is to seek the ensemble mean (u) as discussed in the static context in Section II,B,2 and to hope that this relates in some reasonable, but unspecified way to what would be observed if some ‘local average” were taken. Spatial variation of the mean field is an essential feature of wave propagation and the volume averaging approach given for static problems in Section 11,BJ is simply not available. For applications to composites, for which a transducer is bound to register displacement averaged over a region large relative to the microscale, knowledge of the mean field (u) is likely to be sufficient. For geophysical applications, however, an individual seismogram registers displacement over a scale small relative to the scale of the relevant inhomogeneities and information on fluctuations about the mean is also of interest. A discussion is given by Hudson (1968), for example. Attention will be restricted to the mean wave in this article for the reason that fluctuations have not yet been considered within the very recently developed framework that will be described. Before proceeding with this, brief mention will be made of the longer established methods, which amount to direct adaptations to the elastic case of methods developed for acoustic or electromagnetic waves. First, for a weakly inhomogeneous medium, the displacement field can be described by the appropriate dynamical variant of Eq. (3.6) and this can be solved by perturbation theory. Knopoff and Hudson (1964) and Hudson (1968) performed the most elementary iteration (the Born approximation) to obtain local results that were not valid uniformly over large propagation distances. Karal and Keller (1964) also retained only terms of low order, but obtained

Overall Properties of Composites

65

an integrodifferential equation for (u) by employing the method of smoothing (described in Section 111,AJ) and eliminating uo in favor of (u). Plane-wave solutions were shown to satisfy an eigenvalue problem which gave the wavenumber as a complex-valued function of frequency. Thus, the waves were predicted to be dispersive and also to decay exponentially with distance of propagation. This does not imply physical dissipation since only the mean wave is described in this way: as the wave propagates “coherent” energy in the mean wave is scattered into “incoherent” form by the inhomogeneities. The results of Karal and Keller (1964) are restricted to weak inhomogeneities,but are not limited in frequency. McCoy (1973) considered perturbation series to arbitrary order, obtained by the method of smoothing exactly as in Section III,A,l, to give a formally exact equation governing (u). He then restricted attention to long waves and retained only terms of low order in the ratio microstructural dimension to wavelength. This provided a demonstration that long waves propagate as though the material were uniform, with moduli equal to the overall moduli given by (3.21) and density equal to the mean density of the composite. McCoy also gave the lowest order real and imaginary corrections to the long-wavelength,low-frequency limit of the dispersion relation. An interesting point of principle was settled in this way, but quantitative prediction would still require the summation of series like (3.21). For materials comprising a matrix containing a dispersion of inclusions or fibers, an explicit “multiple scattering” approach is usually adopted, in which the total field is expressed as the sum of an “incident” field uo and fields from each of the scatterers. In calculating the field scattered from inclusion A, the field incident upon A is taken as the sum of uoand the fields scattered from all other inclusions. An infinite set of integral equations (one for each inclusion) is generated by this procedure. It is reduced by taking the expectation value of the equation for inclusion A , conditional upon that inclusion being fixed. This introduces the conditional mean field (u)~, but also, analogously to (4.23), amean field ( u ) ~conditional ~ upon inclusions A and B being fixed. The usual method of closing the heirarchy is to employ Lax’s (1952) quasicrystalline approximation in some form. Datta (1977) took essentially ( u ) ~= ( u ) ~and ~ found that long waves did not propagate as though the material were uniform with overall static moduli and mean density, as predicted by McCoy (1973). However, subsequently (Datta, 1978), he adopted a quasicrystalline approximation on parameters defining the source that generates the field scattered from any chosen inclusion and confirmed McCoy’s finding. Bose and Ma1(1973,1974) simplifiedcalculations by making a “point-scatterer” approximation and again found wave speeds and attenuation coefficients in broad agreement with McCoy (1973). Varadan et al. (1978) and Varadan and Varadan (1979) performed more extensive

66

J. R. Willis

calculations for the propagation of SH waves in a composite containing aligned fibers by expressing both incident and scattered fields as eigenfunction expansions, whose coefficients were related through the T-matrix formalism of Waterman (1969, 1976). The equations were closed by making a quasicrystalline approximation for the coefficients defining the scattered field and solved by truncating the series after several terms. All of this work relies quite heavily upon explicit solutions for point scatterers. It has, consequently, been developed only when the matrix and inclusions are isotropic. The remainder of this section is devoted to an alternative approach, developed by Willis (1980b,c),which comprises a direct generalization of the methods given in Section IV. This does not rely upon isotropy and, although it has so far been applied only to a matrix containing inclusions, it is equally applicable to a composite such as a polycrystal. Thus, it offers a unified approach to wave propagation in any composite. A. POLARIZATION FORMULATION

1. Integral Equations

The equation of motion for any continuum subjected to body force f per unit volume may be given in the form diva + f = p,[,

(8.1)

generalizing (2.17), where p denotes momentum density. Equation (8.1) is independent of material properties. These enter through the constitutive relation (2.16)for the stress and the additional relation

P = PU,t

(8.2)

for the momentum density, where p denotes the mass density of the medium and u is its displacement. Usually, (8.2) is substituted directly into (8.1), but there is some advantage here in treating (8.2) as a constitutive relation, on a par with (2.16).Relative to a comparison material with moduli Lo and density po, the stress polarization z has already been defined by (3.24). The momentum polarization 1 ~ : is now introduced through the corresponding definition (8.3) x = ( P - PO)U,t. Thus, in terms of Lo and p o , a = Loe

+z

and

p = p o ~ ,+t x .

(8.4)

Substitution of these into (8.1) then gives an equation for the displacement in the comparison material, when it is subjected to an extra body force div z -

Overall Properties of Composites

67

z,, . Willis (1980b) showed that the displacement field u so generated could be

expressed in the form, analogous to (3.7),

u = U O - St - Mn,

(8.5) where the field uo is the solution of the boundary value problem for the comparison material in the absence of t and n and S and M are operators related to the Green’s function G for the comparison body: S: zij(x,t ) -+ Jdt’

M: q(x, t )

-+

Jvdx’SPij(x,t, x’, t’)rij(x’,t’),

(8.6)

Jdt’ Jv dx‘Mpi(x,t, x‘,t‘)ni(x’,t‘),

where the kernels of the operators S and M have components SPij(x,t, x’, t’) = dGPi(x,t, x’, t ’ ) / d ~ > I ( ~ ~ ) ,

(8.8)

M J x , t, x’, t’)= - dGpi(x,t, x’, [‘)/at’.

(8.9) It is important to note that the field uo is supposed to satisfy initial conditions of prescribed displacement and prescribed momentum rather than velocity; otherwise, an additional term would appear. Substitution of (8.5) into the definitions (3.24), (8.3) of t, II: now gives the pair of operator equations

+ S,r + M,n = eO, ( p - po)-’n + S,tt + M,+G = ,:u (L - Lo)-’t

(8.10) (8.11)

where e0 is the strain associated with uo and S,, M, have kernels (Sx)pqij= d2Gpi/(dxqdx>)l(ij)(pq), (8.12) (Mx)pqi

= - aZGpi/(dxqW l ( p q ) .

(8.13)

Equations (8.10) and (8.1 1) constitute a generalization of (4.1). Korringa (1972) formulated time-reduced dynamic problems in terms of stress polarization, but did not introduce momentum polarization. Time-reduced versions of (8.10) and (8.11) are obtained by Fourier transformation with respect to time. Now consider a boundary value problem for a random medium which, for definiteness, will be taken as a matrix containing a dispersion of inclusions. As in earlier sections, the matrix will be taken to have moduli L, and density p z . The inclusions, which occupy regions V,, A = 1 , 2 , . . . , have moduli L, and density p1 and are distributed at volume concentration c I . The object will be to determine the mean wave (u) which may be expressed, with the aid of (8.5), in the form (u) = u0 - S(z) - M(a).

(8.14)

68

J . R . Willis

Following the pattern set in Section IV,B, the comparison material is chosen identical with the matrix. The definitions (3.24), (8.3) of the polarizations r, n then show that they differ from zero only over the regions V, occupied by the inclusions and it is useful to define rA,nA as their restrictions to V,. Substitution into (8.10) and (8.11) gives equations similar to (4.19) whose conditional expectations, with V, fixed, may be expressed in forms similar to (4.46):

+

(L, - ~ ~ ) - 1 r : :s,zj

+~

+

, n j JdxB(pB1, - P , ) ( S , ~ ;

+ M,Z;)

+ J’~XBPBIA[S~(~BAB - r;) + M,(~!L- n;)] = eo

+ ~,n;),

- J~X,P,(S,~; (P1

(8.15)

+ S , d + M,PA, + j d x B ( p B I , - P B ) ( s , ~ G + M,~G)

- P2)

-1

A

+ J~x,P,~,[s,~(~:B

-

= upr - J ~ X B P B ( S , ~ ~ ;

+ M,,(~B,,-

+ ~ , ~ n ; ) ,x E v,.

(8.16)

Now it is easy to show that (t) (x, t ) =

j’

d x B p B ~ ( x t), ,

(8.17)

with a similar relation for (n). Thus, from (8.5), the right sides of (8.15) and (8.16) are the mean strain ( e ) and mean velocity (u),,, precisely. Equations (8.15) and (8.16) are not useful as they stand because they involve &, n2B as well as r; , n;. If, however, the quasicrystalline approximation (4.51) is made, for n as well as r, they become approximate but explicit equations for r;, n;. As discussed in Section IV,B,3, the approximation is valid strictly only at low concentrations of inclusions but its alternative derivation in the static case, in Section IV,C, from the Hashin-Shtrikman variational principle, provides some limited encouragement to proceed with it at any concentration. 2. Plane Waves Equations (8.15) and (8.16) apply to any boundary value problem whose solution for the comparison body is uo. The simplest problem in elastodynamics is, however, to seek plane-wave solutions of the differential equations. These waves are usually considered as propagating in an infinite medium. Boundary conditions are necessarily disregarded. The corresponding procedure for Eqs. (8.15) and (8.16) is to set uo = 0 and to relate the operators S, M to the infinite-body Green’s function. Solutions of the homogeneous

69

Overall Properties of Composites

equations are now sought in the form

- x) + wt]},

zi(x,t) = zl(x - x,)exp{-i[k(n ni(x, t) = nl(x - x,)exp{

- i[k(n

x) + ot]},

(8.18) (8.19)

in which the circular frequency w and the unit wave normal n are chosen and the wavenumber k is to be determined. In the same way that elementary plane-wave solutions are defined for a uniform medium, the composite is taken to be statistically uniform. With the forms (8.18) and (8.19), it follows that (z)(x, t ) = clzl exp{ - i[k(n

- x) + ot]}

(8.20)

with a corresponding equation for (n), where TI denotes the mean value of zl(x - xA)over V,. Then, from (8.5) with uo = 0, (u)(x, t ) = -c,[&,

+ MEl] exp{ -i[k(n

*

x) + oil},

(8.21)

where S, M denote the Fourier transforms of S, M, evaluated at (kn,o). Equation (8.21)gives the right sides of the homogeneous equations (8.15) and (8.16), by differentiation with respect to either x or t. The terms involving &, nBABare eliminated by use of the quasicrystalline approximation and the terms that remain are simplified by noting that the time integration in (8.6) or (8.7) produces simply the corresponding Fourier transforms with respect to time; that is, it reduces the operators to time-reduced form. To proceed further, explicit representations of the operators are required. The time-reduced infinite-body Green's function for the comparison material has components Gipthat satisfy the equations (L2)ijklGkp,jl

+ pZwZGip+ 6ip6(x)

=

(8.22)

and it follows immediately by Fourier transforming that

G(kn, w ) = [k2~,(n)- p 2 w 2 ~ ] - ' ,

(8.23)

where L,(n) is the acoustic tensor with components [L2(n)]ik

= (L2)ijklnjnl;

(8.24)

its eigenvalues Ar = p2c? define the speeds c,(n) with which plane waves may propagate in the direction n and the corresponding eigenvectors ur(n) give the wave polarizations. The derivatives with respect to x or t that are required to generate S, S,, etc., from G correspond to multiplication by kn or w for S, S,, etc. A representation for G itself is found by starting from the plane-wave decomposition (4.25)of the delta function. It was shown by Willis(1980b)that

70

J. R. Willis

the eigenvectors u' being taken as orthonormal. Equation (8.25) reduces to (4.29) when w = 0. The low-frequency limit of the dispersion relation for plane waves may now be studied by retaining only terms of order zero in o or k in the homogeneous equations corresponding to (8.15), (8.16). The operators M,, Mtthat appear on the right sides are homogeneous functions of degree zero and so are retained. S , itself, operating on r j or T;, is replaced by the static operator I?", from which it differs by a term of order 0': this follows directly from (8.12) and (8.25). The term of order w 2 in the integral with respect to x, has a bounded coefficient because the integrand contains the factor P,,, - P, which guarantees convergence. The operators M,, S,,, M,, are all at least of order w and so may be disregarded. Thus, to zeroth order, with the quasicrystalline approximation, Eqs. (8.15) and (8.16) reduce to

s,,

s,,

(L,

-~

+ r m r j + Jdx,(~,,, - P,)rmT; = (e),

~)-1rj (p1

- pz)-ln?

= -iw(u),

x

E

v,

(8.26) (8.27)

where (u) is given by (8.21) and (e) has components (eij) = -ikni(uj)/(ijj.

(8.28)

Furthermore, to lowest order, the right-hand sides of (8.26) and (8.27) may be taken constant over V,. Equation (8.26) is identical to (4.52), with z replaced by the local mean strain (e), since the operator r has the same effect as rm on a function such as z; which has bounded support. Given, therefore, that (4.52) generates the estimate f, of the overall moduli, the solution of (8.26) by definition implies (T) =

(L - L,)(e).

(8.29)

Also, if ij is defined as the mean density of the composite, so that p" = ClPl

+ (1 - C J P 2 ,

(8.30)

Eq. (8.27) can be placed in the similar form

(n) = - i o ( p - p&u).

(8.31)

Equations (8.29) and (8.31) show that the mean polarizations (T), (n) are those that would be generated if a wave (u) propagated in homogeneous material with moduli and density p. This agrees with the finding of McCoy (1973).The result can be deduced algebraically by eliminating (2) and (n) between Eqs. (8.21), (8.29),and (8.31).

Overall Properties of Composites

71

Corrections to the low-frequency limit of the dispersion relation are found by retaining some allowance for the exponential in the integrand of (8.25) which defines G. The lowest order imaginary correction is of particular interest, since this determines the rate of decay of themean wave. It is obtained essentially by approximating the exponential that appears in S,, M,t by 1. These constant terms in S,, M,t, when applied to r i , ni, have the effect of producing the meansT,, E l , over V’, premultiplied by the constants and the volume V, of the inclusion. The result of retaining just these terms in addition to those of lowest order in (8.15) and (8.16) is to generate the equations

= (e) (PI

- io3AV, AS,T,,

- pz)-’n;

= -~o(u)

+ ~ c o ~ AAM?rI, V~

(8.32) (8.33)

where ASx, AM are constant tensors defined in the work of Willis (1980b) and the statistics of the inclusion distribution enter the perturbation just through the factor A =1

+ Jdx,(P,la

-p ~ ) .

(8.34)

Equations (8.32) and (8.33) were derived by Willis (1980c), who deduced from them the correction (8.35) (k2/k?)- 1 = (iAM,)Q , , where k, is the zeroth-order approximation to the wave number, n, represents the number density of inclusions and Q, has the form of a scattering cross-section for an inclusion perturbing a wave of rth type. The detailed formulas are given by Willis (1980~) and are not quoted. It may be remarked, however, that the “cross-section” Q, reduces exactly to the scattering crosssection of an isolated inclusion embedded in the matrix, in the limit of low concentrations. The result (8.35) then reduces to the widely accepted form given by Waterman and True11 (1961),so long as A = 1. A similar conclusion has been reached, for acoustic and electromagnetic waves propagating in a medium with isotropic phases, by Twersky (1977, 1978). Equation (8.34) shows that A = 1 - 8 ~ for 1 the well-stirred approximation used, for a corresponding two-dimensional problem, by Varadan, et al. (1978).This approxwhich would correspond to growth imation gives, however, A < 0 for c1 > i, rather than decay of the wave. This could imply failure of the quasicrystalline approximation, but the discussion of Section VI1,C has already shown that the well-stirred approximation is untenable at high concentrations. Twersky (1975) has evaluated A for a Perms-Yevick distribution of spheres : this remains positive at all concentrations.

72

J. R . Willis

The main conclusion of this section is that plane waves propagate as predicted by McCoy (1973) at low frequencies, but, in contrast to the work of McCoy, the approach leads directly to useful estimates of wave speeds and overall moduli that are actually those predicted from the HashinShtrikman formalism. Such estimates are given, together with scattering cross-sections, for composites containing aligned spheroids, by Willis (1980~). It will be apparent, too, that polycrystals could be dealt with in a similar way, following the static reasoning of Section IV,B,3, and self-consistent estimates could be developed by extending the arguments of Section V. Such calculations are in progress at the moment. The problem that remains is to complete the analogy with elastostatics by relating the results to a variational principle. At the time of writing, this has not been accomplished. However, a variational principle, whose implications have not yet been explored, has just been derived; this is outlined in the section that follows. B. VARIATIONAL PRINCIPLE This article is concluded with a derivation of a variational principle, akin to Hamilton’s principle, that is associated with Eqs. (8.15) and (8.16). It is obtained by the same type of reasoning that was applied in Section IV,A to generate the Hashin-Shtrikman principle. First, let b l , e l , p1 be the stress, strain, and momentum fields associated with the displacement u1 generated by polarizations zl, zl,so that ~1

= - S Z ~- M z 1 ,

nl = Loel P1

+ zl,

= POUlJ

+ Zl?

(8.36) (8.37) (8.38)

and define stress, strain, momentum, and displacement fields a2,e 2 ,p2, u2 similarly. It will be assumed, for simplicity,that (8.15)and (8.16)are associated with a boundary value problem with prescribed displacements. Then, u1 = u2 = 0 over aV and the analog of the virtual work equality (4.2) is (8.39) b1,ez) + ( P l , t , U 2 ) = 0, the inner products being defined as spatial means, as in (3.52).It follows now that

J; dt{(z,,ez) -

(7bU2,Jl

Overall Properties of Composites

73

by expressing zl,n1 in terms of cl,p1 as in (8.4) and using (8.39). This is the dynamic analog of (4.5). For any z, n, a Lagrangian functional 9 is defined by the equation

' + S&) +

2 9 = (z,[(L - Lo)--

(7, M,n)

- Po)-' + M J l 4 - (.,S,,t) - 2(2,eo) 2(n,u:). - (n,[(P

+

(8.41)

This generalizes the right side of (4.6). A form of Hamilton's principle for the functional (8.41) may be obtained by considering the variation of the "action" J9ddt. It follows, with the aid of (8.40),that

where Y =

-SZ

Q = POVJ

- Mn

(8.43)

+ a.

(8.44)

The integral in (8.42) vanishes when z, n satisfy (8.15) and (8.16);thus,

6

jd'9 d t = l[ (Q, w 2

- (v,

Wl'd.

(8.45)

Formally, therefore, the requirement that the variation (8.42) should vanish generates the operator Eqs. (8.15) and (8.16).The usual form of Hamilton's principle is expressed in terms of displacement and variations are only allowed which vanish at times 0 and tl . This is reflected by the right side of (8.45) which, however, is less easily reduced to zero by choice of 67, Sn. It is remarked, finally, that the stationary value of the action functional is expressible in the form

ji'9 d t = Ji' (9- go)& - $[(p + po,u - uo)]b',

(8.46)

where 9, are the classical Lagrangian functionals, given by

2 2 = (P,u,,)- (c,e), 2P0= (PO, - (ao,eo),

(8.47) (8.48)

the fields with superscript 0 relating to the solution of the boundary value problem for the comparison material.

74

J . R. Willis

IX. Recent Developments There have been a number of developments, particularly in the area of wave propagation, since this article was completed. For present purposes, the most noteworthy is that a variational structure, related to that indicated in Section VIII,B, is now fully developed (Willis, 198Od). Also, wave propagation has been studied, using the quasi-crystallineapproximation coupled with quantum-mechanical formalism, by Devaney (1980). His work can be related to that described in Section VIII,A, by noting that his “transition operator,” when applied to the field uo in (8.5), generates a source term (div t - z,*).The novelty of Devaney’s scheme is that he proposes to estimate the dynamic Green’s function self-consistently which amounts, in present notation, to defining operators Lo, po, relative to which (t) = (z) = 0. This generalizes the static prescription given by Eq. (5.16). Devaney did not relate his formulation to a variational principle, but this could now be done using the variational principles developed by Willis (1980d). REFERENCES Asaro, R. J., and Barnett, D. M. (1975). The non-uniform transformation strain problem for an anisotropic ellipsoidal inclusion. J. Mech. Phys. Solids 23, 77-83. Batchelor, G. K. (1974). Transport properties of two-phase materials with random structure. Annu. Rev. Fluid Mech. 6, 227-255. Batchelor, G. K., and Green, J. T. (1972). The determination of the bulk stress in a suspension of spherical particles to order c2. J. Fluid Mech. 56, 401-427. Bensoussan, A., Lions, J. L., and Papanicolaou, G. (1978). “Asymptotic Analysis for Periodic Structures,” North-Holland Publ., Amsterdam. Beran, M. J. (1968). “Statistical Continuum Theories,” Wiley (Interscience), New York. Beran, M. J., and McCoy, J. J. (1970). Mean field variations in a statistical sample of heterogeneous linearly elastic solids. Int. J. Solids Struct. 6,1035-1054. Beran, M. J., and Molyneux, J. (1966). Use of classical variational principles to determine bounds for the effective bulk modulus in heterogeneous media. Q. Appl. Math. 24,107-1 18 Bose, S . K., and Mal, A. K. (1973). Longitudinal shear waves in a fiber-reinforced composite. Int. J. Solids Struct. 9, 1075-1085. Bose, S . K., and Mal, A. K. (1974).Elastic waves in a fiber-reinforced composite. J. Mech. Phys. Solids 22, 217-229. Brailsford, A. D. (1976).Diffusion to a random array of identical spherical sinks. J . Nucl. Mafer. 60,257-278. Brinkman, H. C. (1947).A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles Appl. Sci. Res., Sect. A 1, 27-34. Budiansky, B. (1965). On the elastic moduli of some heterogeneous materials. J . Mech. Phys. Solids 13, 223-227. Budiansky, B. (1970). Thermal and thermoelastic properties of isotropic composites. J . Compos. Mater. 4, 286-295. Budiansky, B., and OConnell, R. J. (1976). Elastic moduli of a cracked solid. Int. J . Solids Struct. 12, 81-97.

Overall Properties of Composites

75

Bullough, R., and Hayns, M. R. (1978). Continuum representation of the evolving microstructure prevailing during the irradiation of crystalline materials. In “Continuum Models of Discrete Systems” (J. W. Provan, ed.), pp. 469-502. Univ. of Waterloo Press. Carlson, D. E. (1972). Linear thermoelasticity. In “Encyclopedia of Physics” (C. Truesdell, ed.), Vol. VI a/2, pp. 297-345. Springer-Verlag, Berlin and New York. Chen, H. S., and Acrivos, A. (1978a). The solution of the equations of linear elasticity for an infinite region containing two spherical inclusions. Int. J . Solids Struct. 14, 331-348. Chen, H. S., and Acrivos, A. (1978b). The effective elastic moduli of composite materials containing spherical inclusions at non-dilute concentrations. Int. J. Solids struct. 14,349-364. Chernov, L. A. (1960). “Wave Propagation in a Random Medium” (Engl. transl. by R. A. Silverman). McGraw-Hill, New York. Childress, S. (1972). Viscous flow past a random array of spheres. J. Chem. Phys. 56,2527-2539. Cole, J. D. (1968). “Perturbation Methods in Applied Mathematics.” Ginn (Blaisdell), Boston, Massachusetts. Datta, S. K. (1977). A self-consistent approach to multiple scattering by elastic ellipsoidal inclusions. J . Appl. Mech. 44,657-662. Datta, S. K. (1978). Scattering by a random distribution of inclusions and effective elastic properties. In “Continuum Models of Discrete Systems” (J. W. Provan, ed.), pp. 11 1-127. Univ. of Waterloo Press. Dederichs, P. H., and Zeller, R. (1973). Variational treatment of the elastic constants of disordered materials. Z . Phys. 259, 103-116. Devaney, A. J. (1980). Multiple-scattering theory for discrete, elastic, random media. J. Math. Phys. 21,2603-2611. Eshelby, J. D. (1957). The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R. SOC.London, Ser. A 241, 376-396. Eshelby, J. D. (1961). Elastic inclusions and inhomogeneities. In “Progress in Solid Mechanics” (1. N. Sneddon and R. Hill, eds.), Vol. 11, pp. 87-140. North-Holland Publ., Amsterdam. Fokin, A. G. (1979). Statistical properties of inhomogeneous solid media: Central moment functions of material characteristics. J. Appl. Math. Mech. (Engl. Trunsl.) 42, 570-578. Fokin, A. G., and Shermergor, T. D. (1969).Calculation ofeffective elastic moduli ofcomposite materials with multiphase interactions taken into account. J . Appl. Mech. Tech. Phys. (Engl.Trans(.) 10, 48-54. Gel’fand, I. M., and Shilov, G. E. (1964). “Generalized Functions,” Vol. 1. Academic Press, New York. Gubernatis, J. E., and Krumhansl, J. A. (1975). Macroscopic engineering properties of polycrystalline materials: Elastic properties. J. Appl. Phys. 46, 1875-1883. Hashin, Z., and Shtrikman, S. (1962a). On some variational principles in anisotropic and nonhomogeneous elasticity. J . Mech. Phys. Solids 10, 335-342. Hashin, Z . , and Shtrikman, S. (1962b). A variational approach to the theory of the elastic behaviour of polycrystals. J . Mech. Phys. Solids 10, 343-352. Hashin Z., and Shtrikman, S. (1963). A variational approach to the theory of the elastic behaviour of multiphase materials. J. Mech. Phys. Solids 11,127-140. Hershey, A. V. (1954). The elasticity of an isotropic aggregate of anisotropic cubic crystals. J . Appl. Mech. 21, 236-240. Hill, R. (1952). The elastic behaviour of a crystalline aggregate. Proc. Phys. SOC., London, Sect. A 65, 349-354. Hill, R. (1963a). Elastic properties of reinforced solids: Some theoretical principles. J. Mech. Phys. Solids 11, 357-372. Hill, R. (1963b). New Derivations of some elastic extremum principles. In “Progress in Applied Mechanics,” The Prager Anniversary Volume, pp. 99-106. Macmillan, New York.

16

J . R. Willis

Hill, R. (1964). Theory of mechanical properties of fibre-strengthened materials. I. Elastic behaviour. J . Mech. Phys. Solids 12, 199-212. Hill, R. (1965a). Continuum micromechanics of elastoplastic polycrystals. J . Mech. Phys. Solids 13, 89-101. Hill, R. (1965b). A self-consistent mechanics of composite materials. J . Mech. Phys. Solids 13, 21 3-222. Hinch, E. J. (1977).An averaged-equation approach to particle interactions in a fluid suspension. J . Fluid Mech. 83, 695-720. Hoenig, A. (1979). Elastic moduli of a non-randomly cracked body. Int. J . Solids Struct. 15, 137-154. Hori, M., and Yonezawa, F. (1974). Statistical theory of effective electrical, thermal and magnetic properties of random heterogeneous materials. 111. Perturbation treatment of the effective permittivity in completely random heterogeneous materials. J . Math. Phys. 15, 2177-21 85. Howells, I. D. (1974). Drag due to the motion of a Newtonian fluid through a sparse random array of small fixed rigid objects. 1. Fluid Mech. 64,449-475. Hudson, J. A. (1968).The scattering of elastic waves by granular media. Q. J. Mech. Appl. Math. 21,487-502. Ishimaru, A. (1978). “Wave Propagation and Scattering in Random Media,” Vol. 2. Academic Press, New York. Karal, F. C., and Keller, J. B. (1964). Elastic, electromagnetic and other waves in a random medium. J . Math. Phys. 5 , 537-547. Khatchaturyan, A. G. (1967). Some questions concerning the theory of phase transformations in solids. Sou. Phys.-Solid State (Engl. Transl.) 8, 2163-2168. Kneer, G. (1965). Calculation of elastic moduli of polycrystalline aggregates with texture. Phys. Status Solidi 9, 825-838. Knopoff, L., and Hudson, J. A. (1964).Scattering of elastic waves by small inhomogeneities. J . Acoust. Soc. Am. 36, 338-343. Korringa, J. (1972). Propagating modes of a heterogeneous, macro-homogeneous continuum. J . Phys. (Orsay, Fr.) 33, C6, 117-122. Korringa, I. (1973). Theory of elastic constants of heterogeneous media. J . Math. Ph,ys. 14, 509-513. Kroner, E. (1977). Bounds for effective elastic moduli of disordered materials. J . Mech. Phys. Solids 25, 137 - 155. Laws, N. (1973). On the thermostatics of composite materials. J . Mech. Phys. Solids 21, 9-17. Laws, N., and McLaughlin, R. (1978). Self-consistent estimates for the viscoelastic creep compliances of composite materials. Proc. R. Soc. London, Ser. A 359, 251-273. Laws, N., and McLaughlin, R. (1979). The effect of fibre length on the overall moduli of composite materials. J . Mech. Phys. Solids 27, 1-13. Lax, M. (1952). Multiple scattering of waves. 11. The effective field in dense systems. Phys. Rev. 85, 621-629. Leitman, M. J., and Fisher, G. M. C. (1973).The linear theory of viscoelasticity. In “Encyclopedia of Physics” (C. Truesdell, ed.), Vol. VI a/3, pp. 1-123. Springer-Verlag, Berlin and New York. Lundgren, T. S. (1972). Slow flow through stationary random beds and suspensions of spheres. J . Fluid Mech. 51, 273-299. McCoy, J. J. (1973). On the dynamic response of disordered composites. J . Appl. Mech. 40, 511-517. McCoy, .J. J. (1979). On the calculation of bulk properties of heterogeneous materials. Q.Appl. Math. 36, 137-149.

Overall Properties of Composites

77

Miller, M. N. (1969a). Bounds for effective electrical, thermal and magnetic properties of heterogeneous materials. J . Math. Phys. 10, 1988-2004. Miller, M. N. (1969b). Bounds for effective bulk modulus of heterogeneous materials. J. Math. Phys. 10, 2005-2013. Perkus, J. K., and Yevick, G. J. (1958). Analysis of classical statistical mechanics by means of collective co-ordinates. Phys. Rev. 110, 1-13. Reuss, A. (1929). Calculation of the flow limits of mixed crystals on the basis of the plasticity of mono-crystals. Z . Angew. Math. Mech. 9, 49-58. Rosen, B. W., and Hashin, Z. (1970). Effective thermal expansion coefficients and specific heats of composite materials. Int. J. Eng. Sci. 8, 157-173. Sanchez-Palencia, E. (1974). Comportements local et macroscopique d’un type de milieux physiques hkterogenes. Znt. J. Eng. Sci. 12, 331-351. Schapery, R. A. (1968). Thermal expansion coefficients of composite materials based on energy principles. J . Compos. Mater. 2, 380-404. Talbot, D. R. S., and Willis, J. R. (1980). The effective sink strength of a random array of voids in irradiated material. Proc. R. SOC.London, Ser. A 370, 351-374. Throop, G. J., and Bearman, R. J. (1965). Numerical solutions of the Perkus-Yevick equation for the hard-sphere potential. J. Chem. Phys. 42,2408-241 1. Uscinski, B. J. (1977). “The Elements of Wave Propagation in Random Media.” McGraw-Hill, New York. Twersky, V. (1975). Transparency of pair-correlated, random distributions of small scatterers, with applications to the cornea. J. Opt. Soc. Am. 65, 524-530. Twersky, V. (1977). Coherent scalar field in pair-correlated random distributions of aligned scatterers. J . Math. Phys. 18, 2468-2486. Twersky, V. (1978). Coherent electromagnetic waves in pair-correlated random distributions of aligned scatterers. J. Math. Phys. 19, 215-230. Varadan, V. K., and Varadan, V. V. (1979). Frequency dependence of elastic (SH)-wave velocity and attenuation in anisotropic two-phase media. Wave Motion 1, 53-63. Varadan, V. K. Varadan, V. V., and Pao, Y. H. (1978). Multiple scattering of elastic waves by cylinders of arbitrary cross-section. I. SH waves. J. Acoust. SOC.Am. 63, 1310-1319. Voigt, W. (1889). Ueber die Beziehung zwischen den beiden Elasticitatsconstanten isotroper Korper. Ann. Phys. (Le@zig)[3] 38, 573-587. Walpole, L. J. (1966a). On bounds for the overall elastic moduli of inhomogeneous systems. 1. J . Mech. Phys. Solids 14, 151-162. Walpole, L. J. (1966b). On bounds for the overall elastic moduli of inhomogeneous systems. I f . J. Mech. Phys. Solids 14,289-301. Walpole, L. J. (1969). On the overall elastic moduli of composite materials. J . Mech. Phys. Solids 17, 235-251. Waterman, P. C. (1969). New formulation for acoustic scattering. J. Acoust. SOC.Am. 45, 1417 - 1429. Waterman, P. C. (1976). Matrix theory of elastic wave scattering. J. Acoust. Soc. Am. 60, 56 7 -5so.

Waterman, P. C . , and Truell, R. (1961). Multiple scattering of waves. J . Math. Phys. 2, 512537. Wertheim, M. S. (1963). Exact solution of the Perkus-Yevick integral equation for hard spheres. Phys. Rev. Lett. 10, 321-323. Willis, J. R. (1970). “Asymmetric Problems of Elasticity,” Adams Prize Essay. University of Cambridge, Cambridge, England. Willis, J. R. (1975). The interaction of gas bubbles in an anisotropic elastic solid. J . Mech. Phys. Solids 23. 129-138.

18

J . R . Willis

Willis, J. R. (1977). Bounds and self-consistent estimates for the overall moduli of anisotropic composites. J. Mech. Phys. Solids 25, 185-202. Willis, J. R. (1978). Variational principles and bounds for the overall properties of composites. h “Continuum Models of Discrete Systems” (J. W. Provan, ed.), pp. 185-215. University of Waterloo Press. Willis, J. R. (1980a). Relationships between derivations of the overall properties of composites by perturbation expansions and variational principles. In “Variational Methods in Mechanics of Solids” ( S . Nemat-Nasser, ed.), pp. 59-66. Pergamon, Oxford. Willis, J. R. (1980b). A polarization approach to the scattering of elastic waves. I. Scattering by a single inclusion. J . Mech. Phys. Solids 28, 287-305. Willis, J. R. (1980~).A polarization approach to the scattering of elastic waves. 11. Multiple scattering from inclusions. J . Mech. Phys. Solids 28, 307-327. Willis, J. R. (198Od). Variational principles for dynamic problems for inhomogeneous elastic media. Wave Morion 3, 1- 11. Willis, J. R. and Acton, J. R. (1976). The overall elastic moduli of a dilute suspension of spheres. Quart. J . Mech. Appl. Math. 29, 163-177. Zeller, R. and Dederichs, P. H. (1973).Elastic constants of polycrystals. Phys. Stat. So/idi(b)55, 831-842.