Effective elastic properties for a periodic bicontinuous porous medium

J. Mech. Phys. Solids, Vol. 42, No. 2, PP. 283-305,1994 Copyright 0 1994 ElsevierScienceLtd

Pergamon

Printed in Great Britain. All rights reserved 0022-5096/94%6.00+0.00

EFFECTIVE ELASTIC PROPERTIES FOR A PERIODIC BICONTINUOUS POROUS MEDIUM A.

M. CHAPMAN and

J. J. L. HIGDON

Department of Chemical Engineering, University of Illinois, Urbana, IL 61801, USA (Received 25 January 1993 ; in revisedfbrm

1 September 1993)

AESTRACT THE EFFECTIVE ELASTIC properties of a bicontinuous porous medium are calculated. The model geometry assumes a continuous solid phase consisting of an array of overlapping spheres on a simple cubic lattice. The interstitial spaces form a second continuous phase saturated with a compressible fluid. Following the approach of NUNAN and KELLER (J. Mech. Phys. Solids 32, 259-280, 1984), the anisotropic elasticity tensor C,,, is expressed in terms of three scalar parameters OL,p and y. Numerical computations are performed to evaluate these parameters as a function of volume fraction. The numerical computations are complemented by an asymptotic theory valid for small solid contact area. The variational results of HASHIN and SHTRIKMAN (J. Mech. Phys. Solids 10, 343-352, 1962b) are employed to estimate the elastic properties of isotropic polycrystals formed from randomly oriented grains of the anisotropic medium.

INTRODUCTION EFFECTIVE ELASTIC properties of composite materials are the subject of long standing interest in the mechanics community. In such composites, two or more phases with distinct material properties form a heterogeneous material on the microscopic level, which may be viewed as a homogeneous medium on the macroscopic level. Different types of materials arise depending upon the properties of the individual constituents. A simple composite might be formed from two isotropic linear elastic phases with different bulk and shear moduli. In limiting cases, the medium might consist of an elastic matrix with rigid inclusions or a porous elastic solid in which the moduli of the second phase are zero. On the microscopic level, the individual phases are assumed to be homogeneous, but not necessarily isotropic. For example, polycrystals consist of microscopic grains of varying orientation with each grain having the same non-isotropic elasticity tensor. In the analysis of composite materials, theoretical studies may be divided into three general classes : variational methods, self-consistent approximations and microscopic field calculations. Variational methods apply to the most general class of materials and make no assumptions concerning the microscopic structure of the material. The variational principles yield exact analytical results which provide rigorous bounds for the elastic properties of composites. The major limitation of this approach is that the range between upper and lower bounds may grow large for composites whose constituent phases have markedly different properties. Classic results in this area THE

283

284

A. M.

CHAPMAN and

J. J. L.

HIGDON

include the bounds of VOIGT (19 10) and REUSS (1929) who estimated the effective moduli of polycrystals by assuming that individual grains were in a state of constant strain and of constant stress, respectively. HILL (1952) demonstrated that these estimates provide upper and lower bounds for the true elastic moduli. In effect, the Voigt and Reuss bounds result from a variational approach with the simplest choice of trial functions. HASHIN and SHTRIKMAN (1962a, b, 1963) and HASHIN (1964) describe a general variational method which yields more restrictive bounds by incorporating a broader class of trial functions and extends the variational approach to a wider class of materials. Additional results in this area are presented by WALPOLE (1966a, b, 1969). In self-consistent theories of composites, one considers an inclusion of a given shape embedded in an infinite elastic matrix whose properties equal the effective properties of the medium. The solution of the elasticity equations for the isolated inclusion then yields an implicit equation which may be solved for the effective moduli. The attraction of this approach is that it is readily implemented for a variety of inclusion shapes, often yielding simple algebraic expressions for the effective properties. Predictions from self-consistent theories often show good agreement with experimental data, but such results lack the firm foundation provided by a rigorous solution of the governing equations. For an overview of the use of self-consistent theories for composites, see WALPOLE (1969), CLEARY (1978, 1980) and CLEARY et al. (1980). The third approach to predicting effective properties of composites relies upon the detailed knowledge of the microscope displacement and stress fields. In these theories, one must make specific assumptions concerning the shape of the inclusions and the statistical distribution. A well posed boundary value problem is solved, and the effective properties are evaluated based on integral averages of the microscopic fields. The microscopic field approach provides rigorous results for the effective properties of idealized media. The utility of these results depends upon how well the idealized geometry represents the real systems of interest. The solution of the microscopic field equations is a problem of great complexity and presents a considerable challenge when dealing with large numbers of inclusions. Owing to this difficulty, most efforts have been limited to dilute systems or those with periodic microstructure. RUSSEL and A~RIVOS (1972) analysed elastic materials with dilute suspensions of slender rigid particles. CHEN and ACRIVOS (1978) considered elastic solids with rigid inclusions of spheres including two body interactions to obtain results valid to order c2 where c is the volume fraction of inclusions. For non-dilute systems, NUNAN and KELLER (1984) studied cubic lattices of rigid spheres embedded in an isotropic elastic solid. SANCANI and LU (1987) extended this model to include elastic inclusions of arbitrary moduli, while NEMAT-NASSER and TAYA (1981, 1985), NEMAT-NASSER et al. (1982) and IWAKUMA and NEMAT-NASSER (1983) analysed similar systems and included the effect of non-spherical inclusions. Each of these efforts relied upon numerical computations owing to the complexity of the microscopic displacement fields. In the present effort our interest in the effective properties of composites is associated with the propagation of acoustic waves through fluid saturated porous media. In a classic series of papers, BIOT (1956a, b, 1962a, b) developed a general theory for wave propagation in porous media, addressing such issues as wave speed, attenuation and

Effectiveelastic properties

285

dispersion, and including the effects of elasticity, viscous dissipation and anisotropy. Under the proper set of conditions (CLEARY, 1978 ; BURRIDCE and KELLER, 1981), Biot’s analysis shows that the acoustic response of a medium is characterized by a set of effective properties which may be evaluated independently for the fluid and solid phases. For the fluid phase, the solution of the unsteady Stokes equations determines the dynamic permeability, which is a function of the wave frequency as well as the microstructure of the porous solid. A detailed discussion of this problem is given by CHAPMAN and HIGDON (1992). For the solid phase, the deformation is governed by the elasticity equations for a solid in equilibrium with a fluid at rest, and the effective properties are the elastic moduli of the composite material. Owing to its importance in acoustic wave propagation, the characterization of a fluid saturated elastic solid has received much attention in the composites literature. In addition to the references cited above, significant contributions may be found in the efforts of BIOT and WILLIS (1957), BROWN and KORRINGA (1975), RICE and CLEARY (1976) and CARROLL (1980). These authors discuss the form of the effective elasticity tensors, the relationship between the fluid properties and the effective elastic moduli, and experimental techniques for measuring the effective properties. In the present article, our goal is to calculate the effective elastic properties of a porous solid by solving for the detailed microscopic displacement fields in a well characterized porous microstructure. The most significant feature of the present model is that the fluid and solid phases are interpenetrating and the medium is bicontinuous ; that is, each phase is continuously connected throughout the material. This geometry may be contrasted with that of a suspension, in which one phase is present as a continuous matrix, while the second phase is dispersed in the form of isolated inclusions. A bicontinuous geometry has important consequences for acoustic waves in porous media, because it allows wave propagation in both the fluid and solid phases. The bicontinuous model is representative of the solid microstructure in geophysical media and has relevance for oil exploration, enhanced oil recovery, hydrofracturing, seismic waves and earthquakes. In addition, this morphology is characteristic of synthetic materials such as foams and membranes, which are formed by the solidification of a liquid phase and the subsequent evaporation of a volatile solvent. To generate a simple model of a bicontinuous porous medium, we consider a periodic array of solid spheres on a simple cubic lattice (Fig. 1). As the radius of the spheres is increased beyond the limit of first contact, the overlapping portions are eliminated producing a circular contact area between adjacent truncated spheres.We take this lattice of truncated bodies as our model of a bicontinuous porous medium. This geometry has previously been employed in numerous studies of bicontinuous porous media. SCHWARTZ and KIMMINAU (1987) analysed the electrical conductivity for insulated solids with conducting pores, while CHAPMAN and HICDON (1992) studied unsteady viscous flow through the interstitial spaces. For the elasticity problem, we note that the finite contact area is a characteristic of the bicontinuous porous solid in the undtlformed state. Specifically, we are not considering a Hertzian contact problem where the contact area would depend upon the deformation of the elastic spheres. The present calculations may be considered to be an extension of previous work by Nunan and Keller, Sangani and Lu and Nemat-Nasser et al. All of these studies were restricted to systems with isolated inclusions. Our analysis complements the previous

286

A. M. CHAPMANand J. J. L. HIGDON GRAIN CONSOLIDATION

MODEL

UNIT CELL

FIG. 1. Schematic view of three-dimensional bicontinuous porous medium based on overlapping spheres on a simple cubic lattice. Inset shows unit cell used for computations.

studies and extends the periodic model to new applications involving bicontinuous media. In a companion paper (CHAPMAN and HIGDON, 1992), we have studied the fluid phase and computed the dynamic permeability for this model geometry. With the present results, we calculate the effective elastic properties, and hence provide a complete characterization for the acoustic properties of a bicontinuous porous medium.

FORMULATION

We consider a bicontinuous periodic elastic porous medium saturated with a compressible fluid. The solid phase consists of a system of spheres of radius a centered on a simple cubic lattice with lattice site spacing d and lattice vectors s, = (1 ,O,O)d, s2 = (O,l,O)d, s3 = (O,O,l)d. For radii in the range d/2 < a < d/,,k2, the overlapping spheres are truncated at the cell boundaries to form a continuous solid phase, while the interstitial space forms a second continuous phase saturated with fluid. This configuration constitutes the undeformed equilibrium state of the medium, with the solid contact area determined solely by the model geometry. This truncated sphere model may be contrasted with a Hertzian contact problem in which the contact area between deformed spheres is a function of the applied stress. The volume of the unit

287

Effective elastic properties

cell is VO = d3. The solid phase occupies a region Ds with volume Vs and volume fraction c, while the fluid occupies region DF with volume V, and porosity E E 1 -c. The displacement and stress in the solid phase are u, r, while their counterparts in the fluid phase are Y, O. We assume that the solid is a homogeneous isotropic elastic material with Lame constants I and p. The solid phase strain tensor and stress tensor then satisfy 1 aui

e"=j ( F+z

aZ--

du, >

'

Ti, =

Within the fluid phase, we assume a compressible a fluid at rest, the stress is isotropic and we have cij = -P,,G,~,

2 axj

%,ei,+2Wij,

=

0.

fluid with bulk modulus

p. = -rcF div v.

(~2~3) ICY.For

(4s)

Here (- div v) represents the compression of the fluid when the system is taken from the reference pressure to the current system pressure p,,. We assume that the material is subjected to a homogeneous macroscopic deformation with deformation gradient U,. The displacement is written as the sum of a macroscopic mean component U,x, and a periodic component Ci uj = uiix, + z.&. The macroscopic

(6)

strain tensor is defined as E, = $(U,+

U,).

(7)

In a composite material with two solid phases, it is common to define a mean stress as an integral over the entire volume including both phases. For a fluid filled material, it is useful to define separate averages over the two phases and to define a bulk stress as the sum of these contributions. This definition clarifies the role of the strain Eij and of the fluid pressure p. in contributing to the bulk stress. The mean stress tensors are defined as 1 fij = ~ z,, d V, VS s4

8, = v

1

Oij d V.

F s4

With the constitutive equation (2) and the divergence may be written as a surface integral, while the second leading to 1 [AbijUltI,+ /.L(Uinj+ Ujni)] dS, rii = ~V S s JD, The bulk stress (BIOT, 1962a) is defined as

which -may then be written

&9)

theorem, the first integral may be evaluated directly

~?lj = -PoSi,.

(LO,1 1)

288

A. M. f,i =

Splitting

1

b”0

J. J. L.

1 ~~ ” ~“, [~4,w,+14ujn, V s

U, into its periodic xi, =

CHAPMAN and

HICDON

+ Lp,)] dS-

and non-periodic

Ep&,.

parts leads to

s ?.!I,

(13)

[A6i,ti/n, + ,U(ti,nj + tijni)] dS+ ~(116,,E,I+ ~/_LLE~,) -~pp,6,i.

(14)

To calculate the bulk stress, one might specify the macroscopic strain E,, and the fluid pressure pO. With these specifications, the microscopic displacement field may be computed and the bulk stress evaluated from (14). As a special case, if the solid is “drained”, the fluid pressure p. is zero, and the bulk stress is a function of the deformation E,, alone. In certain applications, the fluid pressure is unspecified, but constraints are placed on the volume of fluid flowing into or out of the medium. For these cases, we follow previous authors (BIOT, 1962a ; BURRIDGE and KELLER, 1981), and introduce < as a measure of the accumulation of fluid within the medium; here, c is based on the displacement of fluid relutiue to the solid. Thus we define (t’i - Ui,.Ylc,)n,dS where the surface of integration is the intersection of the unit cell boundary fluid domain and the normal vector points into the unit cell. The divergence theorem applied to (15) yields [=

-2

o ,,, &(viV s

u,x,,> dV-

k

(u, - U,,xJn, dS. _ 0 s flI,n?D,

(15) with the

(16)

With (5) and (7), the volume integral yields terms E~~/K~ and SE,,. The surface integral covers the fluid-solid boundary where the no penetration condition requires I,‘;= ui ; hence U,- U,x, = ii,. The expression for [ takes the form (17) where the normal points out of the solid. The surface of integration may be taken as the entire boundary of the solid because the integral of z&n,on the solid portion of the cell boundary is zero. Equation (17) is equivalent to (30) in Burridge and Keller. Equations (14) and (17) give the prescription for evaluating the macroscopic quantities ii, and [ once the microscopic displacement field is known. These equations may be viewed as the constitutive equations relating 5 and [ to the imposed macroscopic strain E, and the fluid pressure po. Owing to the linearity of the elasticity equations, these constitutive equations may be written in terms of tensorial coefficients (BIOT, 1962a) as T’i, - = Cllk,Ek,+ 2po

(18)

289

Effectiveelasticproperties

(19) Earlier authors (BIOT, 1962a ; BURRIDCE and KELLER, 1981) have demonstrated the equivalence of the two terms identified as M,/M. These constitutive relationships hold for arbitrary microstructure. For simple cubic lattices, NLJNAN and KELLER (1984) have shown that the elasticity tensor has the form

cl,k, = Cn+ Nlbijskl

+

PL(l+ P)(6k6jl

+

sdsjk>

+

2P(a-B)Si,kl

(20)

where A and ,u are the solid phase Lam& constants and CI,b and y are computed from the microscopic displacement solution. The quantity dilkl = 1 if i =j = k = 1, and is zero otherwise. While Nunan and Keller’s derivation was for an elastic matrix with rigid inclusions, all symmetry arguments carry over to the present circumstances. For cubic lattices of spheres, the remaining tensors Mii are isotropic and we rewrite (IS), (19) as

[ = GE,[+$.

(22)

The parameter G is equivalent to a in the notation of BIOT (1962a), but we have introduced a new symbol to avoid confusion with the parameters associated with C,,,.

PERIODICBOUNDARYVALUE PROBLEMS With the formalism above, we may characterize the effective elastic properties of the medium in terms of five parameters : a, fl, y, G and M. To evaluate these parameters, we solve the governing equations (l)-(3) for the microscopic displacement fields subject to specified macroscopic straining fields. To find the periodic part of the displacement field, we subtract off the macroscopic displacement and its associated stress. First. define uy = U,x, ;

$

= ;lE&j+

2pE,,

(23)

which represent the displacement and stress for a pure solid subject to the imposed strain. Subtracting these quantities from the total displacement and the total stress yields the periodic components 7ii = u,-uup”;

ti, =

zzi-zpp.

To determine the effective properties of the medium, value problems for the periodic displacement Oj.

we consider

(24) three boundary

A. M. CHAPMAN and J. J. L. HIGDON

290

(b)

F~ti. 2. Boundary

value problems for periodic unit cell. (a) Problem 1 : homogeneous strain with E,, = 6, !6,,, (b) problem 2 : homogeneous shear with E,, = 6,, ~3,~ + 6,2d,,

For the first problem, illustrated in Fig. 2(a), let the fluid pressure be zero and the imposed macroscopic strain be a normal strain along the s-axis with E,, = cS,~~,~.Let the surface traction be defined by x = Z,,n,. The boundary conditions on the unit cell are dictated by conditions of periodicity and symmetry. Periodicity requires that the normal displacement be of equal value on opposing faces. On the other hand, symmetry requires equal but opposite values. Combining these results, we conclude that the normal displacement is zero on each face. For the surface traction, symmetry requires that the transverse components be equal on opposing faces, while periodicity requires equal but opposite values. We conclude that the transverse components of surface traction are zero on each face. On the surface of the sphere, the total stress is zero, and we require that the surface traction satisfy ,A = -r; n,. In summary, we have the following conditions for problem 1 : normal strain E,, = 6,,6,, ; pressure p. = 0; and boundary conditions .Yfaces

x = +d/2

6, = 0

.f, = 0

.f = 0

yfaces

I,=$-d/2

.fy=O

ti,. = 0

,K = 0

2 faces

I’ = f d/2

,f; = 0

.r’,. = 0

ii2 = 0

sphere

,f; = -z,T

Y=U

With the specified normal

strain,

”

n,.

(25)

the bulk stress Z”,,in (21) reduces to

f,, = (3.+jQ)fi,,+2~(1

+c()6,,6,,.

(26)

Thus the solution of this boundary value problem combined with (14) yields values for r and y. Next, we consider a boundary value problem with zero fluid pressure, and an imposed shear strain E,, = 6i,6,2 +6r3d,, as illustrated in Fig. 2(b). In this problem, symmetry requires that transverse components of displacement are equal but opposite on opposing faces in the x and 4’ directions. Periodicity requires equal values on opposing faces, and we conclude that the transverse displacement is zero on the x and ?/’faces. For the surface traction, symmetry and periodicity require that the normal component is zero on these faces. On the z faces, the conditions of symmetry and periodicity are as in Problem 1, and we conclude that the normal displacement and

291

Effective elastic properties

transverse surface traction are zero. On the surface of the sphere, the total stress is zero, and the surface traction satisfies A = - rLTnj. In summary for problem 2 : shear strain Ejj = SilSj, +6,,6,, ; pressure p0 = 0 ; and boundary conditions xfaces

x=fd/2

fX=O

zq = 0

ziz = 0

yfaces

y=fd/2

u’,=O

L.=O

ti; = 0

z faces

z = f d/2 fX = 0

_6 = 0

ziz = 0

sphere

j; = -r,Tn,.

Y=U

With the specified shear strain,

(27)

the bulk stress fij in (21) reduces to

crj = 11(1+P)(6ildj2+di28;l)~

(28)

The bulk stress is evaluated using the solution for the microscopic displacement field yielding the value for fi. The third boundary value problem involves zero strain E,, = 0, and an imposed fluid pressure p,, = 1. The solution tii to this problem may be inferred from that of Problem 1 by considering an isotropic compression. Consider an isotropic strain Eli = - iSfjwith fluid pressure p0 = (31”+ 2,~). With these conditions, the displacement satisfies U, E Eijxj exactly and u’, = 0. Comparing expressions (14) and (20), (21) for ?,; yields -(3/?+2/~) Let K = A+?p

= -3(A+/~+j)-2~(1+c~)-G(3A+2~).

be the bulk modulus

of the solid, then (29) reduces to G = -

By similar process,

expressions -3e+

which simplifies

(29)

E(y+$).

(30)

(17) and (22) for i yield

;(3A+2,r)

= -3G+

;(3;+2~)

(31)

to

(32) With results (30) and (32) the parameters G and M are expressed in terms of c( and y together with the physical properties of the individual phases. Given that the evaluation of c( and y relies only on solutions withpo = 0, we conclude that the elastic properties of the fluid saturated media may be inferred from those of the “drained” solid. In the following sections, we present the solution of the microscopic elasticity equations and evaluate the effective elastic parameters.

A. M. CHAPMAN

292

and

J. J. L. HIGU~N

SOLUTION OF ELASTICITY EQUATIONS The solution of the periodic boundary value problems of the previous section might be achieved in a number of ways. The problems could be formulated in terms of integral equations based on the periodic Green’s functions as in the work of Nunan and Keller. Alternatively, one might consider infinite harmonic expansions based on derivatives of the periodic Green’s function as employed by Sangani and Lu. In each of these efforts, great simplification was possible owing to the spherical surface of the inclusions. In the present circumstances, the surface of the inclusion is that of a truncated sphere and these simplifications do not apply. We chose to solve the problem using a collocation method based on Kelvin’s general solution of the elasticity equations. In this approach, the displacement is expanded in a high order series of spherical harmonics, and the boundary conditions are enforced at discrete collocation points on the sphere surface and on the periodic boundaries. By choosing an equal number of unknowns and collocation conditions, one would obtain a large linear system to determine the unknown coefficients. In practice, we have found that such a procedure leads to an ill-conditioned system which is sensitive to the location of the collocation points. To avoid this problem, we employ an excess number of collocation points and solve a least squares system to minimize the errors in the boundary conditions. The governing equations for linear elasticity derive from the combination of the constitutive equation (2) and the force balance (3). The equation for the displacement ii is

(3”+p)v(v-i-l)+pv2ii The general solution of this equation [LOVE (1927), art 1721 in the form ii=

f

in spherical

= 0. coordinates

r’Vo,+~,rw,+V~,,+rxVX,

(33) was given by Kelvin

(34)

where o,, @‘, and xn are solid spherical harmonics. For the geometry of Fig. 1, we choose the origin at the center of the spherical inclusion, Since the origin lies within the solid domain, the solution (34) is restricted to harmonics of positive order with n > 0. Each harmonic is a function of the form x,, = r” i

E(cos

@[A,,,, cos k4 + Dn,k sin k4]

(35)

h = 0

The harmonics o, and Q’, are defined where An,k and D,,, are arbitrary constants. analogously with coefficients (&, En,k) and (c?,,~, Fn,k) respectively. Owing to the symmetry of the spherical geometry and of the boundary conditions, the harmonic expansions for the boundary value problems of the previous section may be substantially simplified. The terms which remain in the expansion are as follows.

293

Effective elastic properties

Problem

Problem

1 : normal

strain

A n,k

for

n=5,7,9,11,...

k=4,8,12,16

,.._

&,k C,,

for

n=0,2,4,6

,...

k = 0,4,8,

D n,k

for

n=3,5,7,9

,...

k=2,6,10,14

,...

En.kF”,k

for

n=2,4,6,8

,...

k=2,6,10,14

,....

12,. .

2 : shear strain

With the symmetry in the spherical harmonics, there is a concomitant reduction in the boundary domain with the variables restricted to the range 4 = 0 to n/4 and 8 = 0 to 7112.In this reduced domain, the boundary surfaces consist of the truncated sphere surface together with a l/S circle on the top plane z = d/2 and a l/4 circle on the front plane x = d/2. Collocation points are distributed over these surfaces with roughly equal density of points on each surface. Details concerning the point distributions are given by CHAPMAN (1993). For each collocation point, the appropriate boundary conditions for displacement or stress are enforced as required by (25) and (27) for the respective boundary value problems. While the collocation method is straightforward in principle, there is one feature which makes the solution more difficult. In the consolidated bicontinuous media of this paper, the surfaces of adjoining spheres intersect at an acute angle. This intersection produces a classic crack geometry with a stress singularity at the line of intersection. The stress singularity leads to slow convergence 0(1/N) in the harmonic expansion, in contrast to the rapid convergence O(e -“) associated with smooth functions. To improve the convergence rate, one might search for simple solutions which match the displacement in the neighborhood of the crack, and thereby subtract from the singularity. Toward this end, consider an axisymmetric body composed of a sphere with top and bottom surfaces truncated at z = *d/2. This is similar to the actual geometry, but without the truncations from the sides of the unit cell. With the boundary conditions (25) for the normal strain problem, the displacement solution for this body is axisymmetric with no 4 dependence. The harmonic expansion for an axisymmetric solution has far fewer terms, and requires proportionately less computation time. Even for the axisymmetric problem, slow convergence will result because the crack is located in the interior of the variable range, i.e. the crack is located at some value 8,,, while 0 ranges from 0 to 71/2. To enhance the convergence, we divided the interval into two subintervals (0, I!?,) and (e,, 7c/2) with separate expansions on each interval. Continuity conditions on stress and displacement are imposed at the border 0 = (3” of the two regions. With this representation, we were able to achieve accurate solutions of the axisymmetric problem employing harmonics up to order N = 150. Even with these high orders, the number of unknowns was quite modest owing to the axisymmetry. Analogous solutions were obtained for truncated axisymmetric bodies corresponding to the x and y periodic cell boundaries. The three high order axisymmetric solutions were appended to the full threedimensional expansion with adjustable weighting coefficients. The ill-behaved displacement field associated with the crack geometry is not axisymmetric, and hence

294

A. M.

CHAIJMAN and

J. J. L.

HEDON

these adjustments will not guarantee a smooth solution. Nonetheless, the convergence of the modified expansion showed dramatic improvement compared with the original expansion. Numerous convergence tests for the axisymmetric solutions and for the three-dimensional expansions are described by CHAPMAN (I 993). As a final point, we note that the stress boundary condition for the normal strain problem ,f = - Z;II, has three non-zero diagonal components (A, 1, I.+ 2~) on the sphere. In the numerical solution, we solved for a single stress component (O,O,l) and built up the required solution by appropriate linear combination. Results for y were also obtained with imposed stress (1.1 ,l) whose solution possesses additional symmetries. For the shear strain, problem 2, there is no simple equivalent of the axisymmetric solutions employed in problem 1. In this problem, therefore, we used the standard three-dimensional harmonic expansion without additional terms. While the basic convergence was limited to 0( i /iv), we found that the errors followed a monotonic pattern allowing extrapolation to higher order results. Examples of this extrapolation are given in the following section.

NUMERICAL RESULTS Using the techniques of the previous section, numerical computations were used to evaluate a, fl and y for a range of volume fractions from 0.53 to 0.90 and for three values of the Poisson ratio v = 0.2,0.3,0.4. The normal strain problem was solved with axisymmetric orders as high as N = 150 and general three-dimensional expansions as high as N = 63. The shear strain problem was solved with expansions as high as N = 91. In all cases, numerous convergence tests were performed as discussed by CHAPMAN (1993). The largest number of unknowns was 783 for the normal strain problem and 1587 for shear strain. In a typical case the number of equations was four times the number of unknowns; however many computations were performed with 10 times the number of unknowns to confirm the insensitivity to collocation point layout. Given the solution for the displacement, Gaussian quadratures were employed to evaluate integrals (14) on the surface of the truncated sphere. Computations were performed on an IBM RS/6000-320 workstation and on a Cray Y-MP. The longest computation times on the Cray ranged from 180 s for the N = 63 normal strain problem to 400 s for the N = 91 shear strain problem. Projections of numerical resolution for selected runs are shown in Table I. The uncertainty in a and ‘4 is inferred from the variation for runs with axisymmetric TABLE I. Numerical accuracy in calculation ofparumeters a, /I and y,for selected concentrationsfor material w’ith Poisson ratio v = 0.2

0.53 0.60 0.80

-0.9448 k 0.0006 -0.7746kO.0010 - 0.4335 + 0.0007

-0.9550+0.007 -0.8114+0.007 -0.4523+0.010

-0.9952&0.001 I -0.9165$-0.0007 -0.531OkO.0026

Effective elastic properties

295

expansions N = 100, 130, 150 and three-dimensional expansions N = 39, 5 1, 63. The predicted uncertainty in fi is based on the difference between the result for N = 91 The and the extrapolated result for N + co. This estimate is extremely conservative. extrapolated values for fi are obtained by fitting values for N = 41, 5 1, 61, 7 1, 8 1, 91 to a linear function in N- ’ . Figure 3 shows the data points and the regression line for three concentrations at v = 0.2. The data for c = 0.60 and c = 0.80 show an excellent fit allowing a clear extrapolation. Based on these data, we believe the actual uncertainty in /I is an order of magnitude smaller than that shown in Table 1. The data for c = 0.53 show more oscillation and only the three highest points (N = 71, 81, 91) were used for the regression line. This concentration is very near the point of close touching spheres (cMIN = 7c/6 z 0.5236), and the lower orders cannot resolve the small contact area between spheres. Even in this case however, the uncertainty is less than 1%. Tables 24 present the values of a, /I and y for a range of volume fractions and for three values of Poisson’s ratio. These parameters represent the difference between the elastic properties of the bulk and the elastic properties of the pure solid. To show the behavior of the bulk properties themselves, the quantities (I+ a), (1 +/I) and (1 + py/A) are plotted in Figs 4-6. These figures show the smooth behavior of the numerical results as they approach the analytical limits at cMIN = 7c/6 and c = 1.0. For further

1’05 5 1.04 _

1.03 _ d c1 1.02 _

0.0

0.005

0.010

0.015

0.020

0.025

l/N FIG. 3. Convergence

linear regression

of parameter b as a function of inverse order of harmonic expansion. Solid lines are curves. Poisson ratio v = 0.2. volume fractions c = 0.53 0, 0.60 0, 0.80 0.

296

A. M. CHAPMANand J. J. L. HEDON TABLE 2. Parameter fraction ,for dtflerent v. Additional entries

dictions

based on asymptotic

V = 0.2 0.5236 0.53 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90

u as function qf volume values of the Poisson ratio marked with an * are pre-

- I .oooo - 0.9448 - 0.9436* - 0.8793 -0.8854* -0.7746 - 0.6883 -0.6078 - 0.5227 -0.4335 -0.3419 -0.2399

v = 0.3

- 1.oooo -0.9366 -0.9335* - 0.8634 - 0.8690* - 0.7500 -0.6617 -0.5751 -0.4917 - 0.4042 -0.3132 -0.2204

theory v = 0.4

- 1.oooo - 0.9262 - 0.9248* -0.8443 - 0.8472* -0.7223 - 0.6280 -0.5423 -0.4577 - 0.3723 -0.2841 - 0.2005

confirmation of the results, we developed an asymptotic theory (see Appendix) valid for small contact areas at c - cMIN c 1. The asymptotic predictions are shown with an * in Tables 24. In all cases, the numerical results smoothly approached the limiting values from the asymptotic theory. While tabulations provide a useful reference standard for the numerical results, algebraic representations are often more convenient for applications. For each of the parameters CC,p and y, we know the exact values at c = cMIN and c = 1.O and the leading order correction for c-cMIN << I. Simple empirical equations which satisfy these constraints and correlate the numerical data take the form TABLE

3. Parameter

,fraction

,for d@erent

fi as jimction values

of volume

of the Poisson

ratio

V. Additional entries marked with an * are predictions based on asJ>mptotic theory 4 ___0.5236 0.53 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90

1’= 0.2

1’= 0.3

L’= 0.4

- 1.oooo -0.9550 - 0.949s* -0.8959 -0.8981* -0.8114 -0.7312 -0.6456 -0.5536 -0.4523 ~ 0.3434 - 0.2268

- I .oooo

- 1 .oooo

-0.9473 -0.9469* -0.8955 - 0.8922* -0.8049 -0.7226 - 0.6360 -0.5447 -0.4439 -0.3364 -0.2222

- 0.9512 -0.9436* - 0.8954 - 0.8854* - 0.8000 -0.7154 -0.6288 -0.5353 - 0.4358 -0.3281 -0.2157

Effective elastic properties

297

TABLE 4. Parameter y as function of volume fraction for d@erent values of the Poisson ratio v. Additional entries marked with an * are predictions based on asymptotic theory

4

v = 0.2

v = 0.3

- 1.oooo

- 1 .oooo -0.9952 - 0.9956* - 0.9772 -0.9818* -0.9165 -0.8387 -0.7408 - 0.6457 -0.5310 -0.3895 - 0.2524

0.5236 0.53 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90

v = 0.4

- 0.9968 - 0.9972* - 0.9852 - 0.9886* - 0.9450 -0.8856 -0.8219 - 0.7376 -0.6358 -0.5067 -0.3532

- 1.oooo -0.9982 - 0.9985* -0.9919 - 0.9938* - 0.9697 -0.9379 - 0.8952 -0.8385 -0.7620 -0.6505 -0.4857

0.6

0.0

1

I

I

I

I

II

0.6

0.7

0.8

0.9

a

1.0

C FIG. 4. Effective

tations,

modulus 1 +a as a function solid line is (36). From bottom,

of volume fraction c. Open symbols are numerical compucurves for Poisson ratios v = 0.2, \I = 0.3 and Y = 0.4.

A. M. CHAPMAN

298

and J. J. L. HIGD~N

0.8

0.6

0.2

0.0 II

0.6

B

0.7

0.8

0.9

1.0

C FIG. 5. Effective modulus I +fi as a function tations, solid line is (37). From bottom.

of volume fraction c’. Open symbols are numerical compucurves for Poisson ratios v = 0.2, v = 0.3 and 13= 0.4.

1 +z = E”‘[U,j(l -?)+(‘+(I 1+B= I+?

‘:“*[h”(l-?)+?+(l

-c’)(a,?+u*?‘)]

(36)

-?)(h,?+h,P)]

(37)

= ;.[y”(l-P)+E+(l--)(y,~+.92S2)]

(38)

where ? = (~-c~,~)/(l -cMIN), and the coefficients LI, b and y are tabulated in Table 5. Note that a”, b, and go are exact values from the asymptotic theory, while the other coefficients are empirical values. The curves represented by (36))(38) are plotted as the solid curves in Figs 4-6. These figures demonstrate that the computed values are well represented by the simple algebraic relationships. Tables 24 and (36))(38) constitute the principal results of this paper. These data provide a complete description of the elastic properties for our model of a bicontinuous porous medium. The parameters a, /I and y characterize a drained solid, while the additional parameters G and A4 (30), (32) are required to describe a general fluid saturated medium. For the special case of an undrained solid with no fluid transfer, we may set i = 0 in (21), solve for p0 in (22) and evaluate the elastic parameters for the undrained medium. The parameters c( and /I’are unchanged, while y satisfies

Effective elastic properties

299

0.8

FIG. 6. Effective computations,

modulus 1+pi/i as a function of volume fraction c. Open symbols are numerical solid line is (38). From top, curves for Poisson ratios 18= 0.2, 1’ = 0.3 and Y = 0.4.

+ G‘M CL+ w)““drai”cd = (2+ whramec~ which is equivalent

(39)

to (3.6) in BIOT (1962b).

EFFECTIVE PROPERTIES OF POLYCRYSTALS Given that the elasticity tensor of our model porous medium is non-isotropic, it is interesting to consider the effective properties of an isotropic polycrystal composed of randomly oriented grains whose individual properties are given by our model’s elasticity tensor C,,,,. Let ,u, = ~(1 +a) and pclg= ,u(I +/I) designate the extensional and shear moduli of the non-isotropic medium and let ti elr = (i*+W)+

Ml

+cI)

(40)

designate its effective bulk modulus. For crystals of cubic symmetry, the bulk modulus of the polycrystal is identical to that of the individual grains; hence ~~~ = KeF. For the shear modulus of the polycrystal, let pLuand pL denote the upper and lower bounds given by HASHIN and SHTRIKMAN (1962b). In the present notation, their results (3.30) and (3.31) may be written

A. M. CHAPMAN and J. J. L. HIGDON

300

TABLE 5. Co@cients,for 1’ =

UC UI a2 h, b, h* g0 9I S? '710 tn, rn2

0.2

0.486767 -0.037885 -0.235127 0.432682 -0.439338 0.461321 0.327951 0.391159 0.105357 0.453510 -0.284270 0.192512

(36)-(38),(44)

v = 0.3

1'= 0.4

0.556306 -0.029772 -0.177130 0.458134 -0.457044 0.490207 0.205274 -0.051243 -0.108973 0.495010 -0.294699 0.240361

0.649023 -0.079439 -0.075239 0.486767 -0.525152 0.597142 0.111070 -0.414492 -0.772510 0.545870 -0.357695 0.354789

(41)

(42) where (43) The numerical results for these bounds are presented in Table 6 as a function of volume fraction and Poisson ratio. Because pX and pL8are of comparable magnitude. the upper and lower bounds are nearly identical, varying only in the fourth decimal place. An empirical equation may be used to correlate these data yielding a prediction for the effective shear modulus of the polycrystal

where the coefficients mo, nz, and nz2 are given in Table 5. The effective properties lcpc and ppc of the isotropic polycrystal are given by equations (40) and (44) with values for cxand y given by (36) and (38). An explicit expression for the Lame parameter i pc may be inferred from the relation lcpc = ~,,r. The result is * = IbPC

(A+w)+:(/A-I*Pc).

(45)

The precise numerical results presented in this paper are based on the idealized geometry assumed in our model. It is interesting to compare the results for the bulk and shear moduli lcpc and ppc with the bounds for an isotropic medium when no restrictions are placed on the microstructure. These bounds have been given by HASHIN and SHTRIKMAN (1963). (4.1)-(4.4). The lower bounds in each case are identically

Effective elastic properties

301

TABLE 6. Efectiue mod& and bounds for isotropic polycrystals v = 0.2 4

1 = 0.6667 G.TIK

K = 1.3333 KMAX /K

ILIP

KJlP

PLIP

&r/P

PMAXIP

0.5236 0.53 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90

0.0000 0.0300 0.0718 0.1545 0.2365 0.3257 0.4158 0.5177 0.6343 0.7538

0.3547 0.3605 0.3793 0.4286 0.4815 0.5385 0.6000 0.6667 0.7391 0.8182

0.0000 0.0552 0.1207 0.2254 0.3117 0.3922 0.4773 0.5665 0.6581 0.7601

0.0000 0.0488 0.1105 0.2026 0.2852 0.3691 0.4585 0.5552 0.6572 0.7680

0.0000 0.0488 0.1104 0.2025 0.2852 0.3691 0.4585 0.5552 0.6572 0.7680

0.0000 0.0450 0.1041 0.1886 0.2688 0.3544 0.4464 0.5477 0.6566 0.7732

0.3547 0.3605 0.3793 0.4286 0.4815 0.5385 0.6000 0.6667 0.7391 0.8182

1’= 0.3 9

i = 1.5000 ‘&R./K

P”IP

/k/P

&/P

hV.4xlP

0.5236 0.53 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90

0.0000 0.0217 0.0523 0.1150 0.1833 0.2540 0.3381 0.4355 0.5528 0.6877

0.2951 0.3005 0.3177 0.3636 0.4143 0.4706 0.5333 0.6038 0.6834 0.7742

0.0000 0.0567 0.1164 0.2155 0.3004 0.3873 0.4759 0.5717 0.6728 0.7785

0.0000 0.0567 0.1162 0.2153 0.3003 0.3872 0.4758 0.5717 0.6728 0.7785

0.0000 0.0527 0.1045 0.1951 0.2774 0.3640 0.4553 0.5561 0.6636 0.7778

0.3654 0.3713 0.3903 0.4400 0.493 1 0.5500 0.611 I 0.6769 0.7480 0.8250

!J = 0.4

i = 4.0000 KC&

li = 4.6667

4

PUIP

PLlP

&/P

0.5236 0.53 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90

0.0000 0.0121 0.0292 0.0656 0.1064 0.1552 0.2159 0.2937 0.40 18 0.5550

0.1963 0.2004 0.2136 0.2500 0.292 1 0.3415 0.4000 0.4706 0.5574 0.6667

0.0000 0.0534 0.1228 0.2282 0.3169 0.4038 0.4944 0.5889 0.6892 0.7904

0.0000 0.0530 0.1225 0.2279 0.3166 0.4036 0.4943 0.5888 0.6892 0.7904

0.0428 0.1046 0.2000 0.2845 0.3712 0.4647 0.5642 0.6719 0.7843

2.1667 KMAX/K

K =

KMAx

I K

0.0000 0.0634 0.1366 0.2500 0.3383 0.4249 0.5083 0.5958 0.6868 0.7796

PJP 0.0000 0.0738 0.1557 0.2777 0.3720 0.4577 0.5423 0.6277 0.7159 0.7995

0.0000

PMAXIP

0.3791 0.3852 0.4044 0.4545 0.5078 0.5645 0.6250 0.6897 0.7589 0.8333

zero, corresponding to any medium in which the solid forms a dispersed phase. The values of the upper bounds ~~~~ and ,uMAXare presented in Table 6. In all cases, we find that the upper bounds are well above the predicted values for the specific microstructure. This confirms our earlier assessment that variational principles are of limited value when the individual phases have markedly different properties.

302

A.

M.

CHAPMAN

and

J. J. L. H~c;oou

ACKNOWLEIXEMENTS This work was supported by the Petroleum Research Fund administered by the American Chemical Society. Additional funding was provided by the National Science Foundation. Computations were performed on the Cray Y-MP at the National Center for Supcrcomputmg Applications at the University of Illinois.

REFERENCES BIOT. M. A. (19.56a) General solutions of the equations of elasticity and consolidation for a porous material. J. Apple. Mech. 23, 91-95. BIOT, M. A. (1956b) The theory of propagation of elastic waves in fluid saturated porous solids, I-&low frequency range, II-high frequency range. J. Acoust. Sot. Am. 28, 168 17X and 179-191. RIOT. M. A. (1962a) Generalized theory of acoustic propagation in porous dissipative media. J. Arwust. Sec. Am. 34, 1256-l 264. Bror. M. A. (1962b) Mechanics of deformation and acoustic propagation in porous media. J. Appl. Php. 33, 1482~1498. BIOT, M. A. and WILLIS, D. G. (1957) The elastic cocflicicnts of the theory of consolidation. J. Appl. Mech. 24, 594-601 BROWN. R. J. S. and KORRINGA. J. (1975) On the dependence of the elastic propcrtics of a porous rock on the compressibility of a pore fluid. Gcop/~,$c.s 40, 608 ~616. BURRIDGE, R. and KELLER, J. B. (1981) Poroclasticity cquatlons derived from microstructure. J. Acoust. Sot. Am. 70, 1140-I 146. CARROLL, M. M. (1980) Mechanical rcsponsc of fluid-saturated porous materials. T/7rowrid and App/ic~tlMec,hrmic,s (cd. F. P. J. RIMROTTand B. TAIIARROE;).North-Holland, Amstcr-dam. CHAPMAX, A. M. (1993) Ph.D. dissertation, University of Illinois. CHAPMAN. A. M. and HIGIION, J. J. L. (1992) Oscillatory Stokes How in periodic porous media. P/I~s. Fluitls A4, 2099- 2 I 16. CHES, H.-S. and A~RIVOS. A. (1978) The effective elastic moduli of composite materials containing spherical inclusions at non-dilute concentrations. Int. J. Solids Struct. 14, 349 364. CLTARY, M. P. (1978) Elastic and dynamic response regimes of fluid-impregnated solids with diverse microstructurcs. ht. J. Sokls Sttwt. 14, 795%? 19. CLEARY, M. P. (1980) Wave propagation in fluid-infiltrated porous media-some rcvicw and analysis. Reports of’ R~~secwch in Mrchanics cd Muterids. MIT, MA. CLEARY, M. P., CHI:K, I.-W. and LEE, S.-M. (1980) Self consistent techniques for heterogeneous media. J. Engng Md., ASCE 106, X61- 887. HASHIN, Z. and SHTRIKMAN, S. (1962a) On some variational principles in anisotropic and nonhomogeneous elasticity. J. Mrch. Ph,vs. Solids IO, 335-342. HASHIN, Z. and SHTRIKMAN. S. (1962b) A variational approach to the theory of the elastic behavior of polycrystals. J. Mech. Phys. Solids 10, 343 352. HASHIN. Z. and SHTRIKMAN, S. (1963) A variational approach to the theory of the elastic behavior of multiphasc materials. J. Mcch. Ph~~s. Solid.s 11, 127-140. HASHIN, Z. (1964) Theory of mechanical behavior of heterogeneous media. .4p/)l. .IMcc/7.Rw. 17, I-9. HILL.. R. (1952) The elastic behaviour of a crystalline aggregate. Phj7.s. SW. Land. A65, 349.m 354. IWAKUMA, T. and NEMAT-NASSER, S. (1983) Composites with periodic microstructure. Cornputrrs Stvuct. 16, 13-20. LOVE, A. E. H. (1927) A Treatise on the Mathematicd Theor), of’ Elasticity. Cambridge University Press. Cambridge. MINDLIN, R. D. (1949) Compliance of elastic bodies in contact. .I. Appl. Mech. 16, 259-268.

303

Effective elastic properties

NEMAT-NASSER,S. and TAYA, M. (198 1) On effective elastic moduli of an elastic body containing

periodically distributed voids. Q. Appl. Math. 39,43-59. NEMAT-NASSER, S. and TAYA, M. (1985) Body containing periodically distributed voids, comments and corrections. Q. Appl. Math. 43, 187--l 88. NEMAT-NASSER,S.. IWAKUMA, S. and HEJAZI, M. (1982) On composites with periodic structure. Mech. Mater. 1,239-267. NUNAN, K. C. and KELLER, J. B. (1984) Effective elasticity tensor of a periodic composite. J. Mech. PhJ,s. Solids 32, 259-280. REUSS, A. (1929) Berechnung der Fliessgrenze von Mischkristallen auf Grund Plastizitgtsbedingung fiir Einkristalle. ZAMM 9,49. RICE, J. R. and CLEARY, M. P. (1976) Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents. Rea. Geophys. Space Phps. 14, 227241. RUSSEL,W. B. and A~RIVOS, A. (1972) On the effective moduli of composite materials : slender rigid inclusions at dilute concentrations. ZAMP 23,434464. SANGANI, A. S. and Lu, W. (1987) Elastic coefficients of composites containing spherical inclusions in a periodic array. J. Mech. Phys. Solids 35, l-2 1. SCHWARTZ, L. M. and KIMMINAU, S. (1987) Analysis of electrical conduction in the grain consolidation model. Geophysics 52, 1402-I 411. SHANKS, D. (1955) Non-linear transformations of divergent and slowly convergent sequences. J. Math. Phys. 34, l-42. TIMOSHENKO,S. P. and G~OIIIER, J. N. (1970) Theory ?fEElasticity. McGraw-Hill, New York. Vorc~, W. ( 19 10) Lehrbuch der Kristallphysik. Teubner, Berlin. WALPOLE, L. J. (I 966a) On bounds for the overall elastic moduli of inhomogeneous systemsI. J. M&. Phys. Solids 14, 151-162. WALPOL~, L. J. (1966b) On bounds for the overall elastic moduli of inhomogeneous systems11. J. Me& Phys. Solids 14, 289-301. WALPOLE, L. J. (1969) On the overall elastic moduli of composite materials. J. Mech. Phys. Solids 17, 235-25 I.

APPENDIX In this Appendix, we develop simple expressions for the parameters c(, fi and y for the asymptotic limit when the spheres just overlap. More specifically, let 6 = c--c~,~ and consider the limit 6 << 1. The radius of the spheres a and the lattice length dare related by d = 2a+O(S). The contact area between adjoining spheres is a circle of radius b where to leading order

Consider

an imposed

normal

stress T” such that the mean stress and macroscopic

f,, = i’i

With these expressions,

0

the constitutive

Z;

;;”

equation

E,,

(20) is evaluated

0 = (~.+/Q~)E,,+~,u(I

To = (~+w)&+a41

Ej.),

strain are

(A.], A.2)

for ?, ,, Z3j yielding

+a)E,,

+“)E,,.

(A.3, A.4)

Solving for CIand y yields (A.5, A.6)

304

A. M.

CHAPMAN and J. J. L.

HIGDON

To find x and ;‘, we need to find the normal strains for a given stress 7”. From a force balance on the solid combined with the definition of mean stress (S), the total force exerted on the top surface of the body is F; = d27,, z 4~~7,. For small S, this force corresponds to the uniform displacement of a circle of radius h normal to a surface which is to leading order an infinite plane. Thus the displacement may be inferred from the classic solution for a rigid circular die impressed in a plane (with zero transverse stress). Note that this is not the Hertz solution for contact between two spheres. The consolidated medium has a finite contact area in the absence of stress, while Hertzian contact gives a contact area which is a function of the applied force. From the solution for the circular die (TIMOSHENKO and GOODIER, 1970. $139), the displacement U, at the top contact surface is 1-V - F; uZ = 4Llh To leading order, the strain E?, = K/U. Substituting

(A.7) for u; and FZ yields

In the absence of horizontal normal stresses, the displacement u, at the front contact surface (I = d//2) arises from the normal stress at the poles. To leading order, this equals the displacement at the equator (0 = n/2) for a sphere subject to point forces F, applied at the poles (0 = 0, 7~). From the l/r decay of a point force, we infer that the displacement U, is smaller than U, by a factor O(hja) and hence E, , c. E,,. Given this result, we combine (A.5) and (AX) and substitute in terms of 6 to obtain (A.9) To evaluate 7, we require a more detailed analysis for u,. LOVE (1927, art. 185) gives the solution for the displacement on a sphere subject to specified surface tractions. Let ,f’= (j;. ,/;. ,J) be the surface traction vector. For point forces k F, applied at the poles, the surface traction becomes (0, 0. .f) where ,J is expressed in axisymmetric surface harmonics as (A.10)

The coefficients d,,, arise from the harmonic expansion for a Dirac delta function on the surface of a sphere. With these coefficients, we evaluate Love’s solution at (0 = n/2, C/J= 0) to find the displaccmcnt u,

(A.11) In this expression,

P,*, (0) is an even order Legendre

polynomial

which is given by

For convenience, let the bracketed term in (A.1 1) be designated S, where the subscript indicates that the term is a function of the Poisson ratio V. This term requires the summation of an infinite series which converges quite slowly : the Nth partial sum has an error order 1IN. Fortunately, this type of oscillating series responds quite well to repetitive application of

Effective elastic properties

305

SHANKS (I 955) transformation. With this approach, 15 terms are sufficient to compute S,, to seven significant figures. The values for v = 0.2, 0.3 and 0.4 are S, = 0.18 11807,0.1953562 and 0.2070935, respectively. The displacement becomes u,

F.

=-s,..

(A. 12)

va

The stram E, , = u,/a takes the form (A.13) Combining

(A.6), (A.8) and (A. 13) yields the expression (A.14)

Substituting

the result (A.9) for (I + 2) and expressing

p/i in terms of v leads to (A.15)

To calculate the parameter and macroscopic strain are

The constitutive

equation

p, consider

an imposed

shear stress such that the mean stress

(20) gives 70 = 2P(l +D)E,,.

(A.18)

The mean stress tensor may be associated with a force couple +F, applied to the faces x = &d/2 and a couple f F, applied to the faces y = &d/2. From the definition of mean stress, the magnitude of the forces is F, = Fv = d%,, z 4~~7,. In the limit as 6 + 0, the force corresponds to the uniform displacement of a circular region of radius b parallel to the surface of an infinite plane wall. MINDLIN (1949) has solved this problem with the result 2-v u, = --F,. e/3

(A.19)

The shear strain is E, z = u,/a which gives

E,2:2-Va3!, 2 Combining

b,u

(A. 18) and (A.20) gives the expression

(A.20)

for b I.‘2

l+8=2+; This completes

the specification

0 3

of the parameters

+0(h). in the limit of small contact

(A.21) area 6 << 1.