The effective mechanical properties of random porous media

J. Mech. Phw. Solids, Vol. 44, No. IO, pp. lYS1620.

Pergamon PI1 SOO22-5096(96)00051-8

THE EFFECTIVE MECHANICAL RANDOM POROUS

1996 Copyright 8 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0022-5096/96 $15.00+0.00

PROPERTIES MEDIA

OF

J. POUTET,tf D. MANZONI,_F F. HAGE-CHEHADE,? C. J. JACQUIN,? M. J. BOUT&A,? J.-F. THOVERTS and P. M. ADLER?8 t lnstitut Fraqais Chasseneuil.

du Petrole, 92506 Rueil-Malmaison Cedex France, 1 LPTMjCNRS, 86 360France, $ IPGP, tour 24,4 Place Jussieu, 75252 Paris Cedex 05. France

(Received 17 October 1995; in reaisedform

25 March 1996)

ABSTRACT The effective mechanical properties of homogeneous porous media can be theoretically determined by a multiple scale expansion. In such a scheme, the media are modelled as being spatially periodic. The elastic equations with appropriate boundary conditions are numerically solved on the unit cell to determine the macroscopic mechanical properties of the medium. Three types of unit cell were investigated : deterministic, fractal and random ; the last type includes media derived from site percolation, random packings of spheres and a real sample of Fontainebleau sandstone. A systematic comparison is made with existing numerical data and with the predictions of various types of renormalization models. Copyright $3 1996 Elsevier Science Ltd Keywords: B. constitutive C. numerical algorithms.

behaviour,

B. granular

1.

material,

B. porous

material,

C. finite differences.

INTRODUCTION

The determination of the deformation under stress of porous materials is a longstanding problem of great interest, for instance, for the oil industry (Bout&a, 1992). The macroscopic stress-strain relationship may be derived either from thermodynamics (Coussy, 199 1) or by homogenization (Sanchez-Palencia, 1980). This upscaling may have two purposes ; as previously mentioned, the first one is to derive the equivalent macroscopic mechanical properties of the rock, which may be useful for large-scale modelling ; the second one is to obtain the local deformations in the pore structure and thus subsequent changes in some important transport properties. such as permeability. In this paper, only the former point is addressed. Various methods have been proposed to derive the macroscopic elastic properties of composite materials such as rocks (see Watt et al., 1976). Only the general references are given in this introduction and more specific references will be provided when necessary. Among the most popular are the variational methods which were initiated long ago by Hashin and Shtrikman (1962a, b). Another classical path is the use of self-consistent approximations which were summarized by Cleary et al. (1980). The last path is the most ambitious since it consists of solving the local elastic equations ; of course, the number of configurations that can be analytically solved is deceptively 1587

1588

J. POUTET ef al.

small and their practical use is limited such as the analysis of a dilute suspension of particles by Russel and Acrivos (1972) and of a cubic array of spheres of Phan-Thien (1983). Hence, one has to use numerical tools to solve general configurations, or at best mixed tools which combine analytical preliminary calculations with numerical computations ; a recent example of such a mixed technique is the paper by Chapman and Higdon (1994) who solved a three dimensional array of overlapping spheres by means of a collocation method based on Kelvin’s general solution of the elasticity equations in terms of spherical harmonics. Finally, let us mention several contributions very closely related to the present one ; Day et al. (1992) and Snyder et al. (1992) determined the elastic moduli of a sheet containing circular inclusions and holes by a discretized-spring scheme ; more recently, Garboczi and Day (1995) used a finite element algorithm to analyse the macroscopic properties of three-dimensional composites with equal phase Poisson ratio. The main purpose of this paper is to provide numerical results for the major possible local structures of porous media. Three such local structures can be distinguished, namely deterministic, fractal and random. They will be successively analysed in this paper. This paper is organized as follows. The general background on porous media is provided by Section 2. The schematization of homogeneous media by spatially periodic structures composed by identical unit cells is presented. The four major local structures, i.e. the contents of the unit cells, are described and illustrated; the generation process is also given, as well as the most important geometrical properties such as the percolation threshold. References to more detailed presentations are made. In Section 3, the elasticity equations are presented. The multiple scale technique is briefly recalled and the effective mechanical properties as well as their major characteristics are introduced. The major features of the numerical method and of the program are described. This program is very versatile in the sense that it was designed in such a way that it is able to handle any geometry in two and three dimensions that can be discretized on a rectangular grid. The resulting system is solved by a classical conjugate gradient method. Various checks were performed on the program, such as symmetry and periodicity checks, which proved to be very efficient for the elimination of programming errors. The influence of the computation precision upon the effective properties of a random medium and of a Menger sponge at generation 2 is analysed. A good quantitative comparison with the semi-analytical results of Sangani and Lu (1987) and Iwakuma and Nemat-Nasser (1983) for cubic arrays of spherical void inclusions is obtained. Section 4 is devoted to deterministic structures such as regular arrays and fractals. Equivalent Young’s modulus and Poisson’s ratios are introduced for materials with one or two symmetry planes; the solids under consideration may be either two- or three-dimensional structures. Power laws, for all the effective coefficients are almost immediately obtained within a few percent. The ratio between the coefficients relative to two successive generations can be derived by various types of renormalization arguments ; for instance, a Sierpinski carpet at generation N can be viewed as being made of carpets at generation N- 1 which behave independently. The predictions of these arguments are compared to the complete numerical results and discussed. The agreement is generally satisfactory, but its quality depends on a crucial way upon the

Random porous media

1589

value of the Poisson ratio. A recent theoretical prediction, which states that the Young modulus of a two-dimensional material containing holes of a given shape is independent of the Poisson ratio of the matrix material (cf. Thorpe and Jasiuk, 1992) is well verified. These computations were extended to the classical three-dimensional structures, which are the fractal foams and the Menger sponges. Computations were also made up to generation 3. Power laws for the effective coefficients are obtained in both cases, but the asymptotic regime is reached almost immediately for fractal foams, in contrast with the Menger sponges where this behaviour is reached for higher generations. Three-dimensional extensions of the previous renormalization arguments are applied ; generally speaking, they are in better agreement with numerical results for fractal foams than for Menger sponges; this might be due to the fact that the former are more isotropic in character than the latter. Various types of random media are addressed in Section 5, which begins with random structures derived from site percolation. The isotropic character of the macroscopic mechanical properties is verified first. Then, systematic computations were performed for various porosities and various sizes of the unit cell. Close to the critical porosity for which the solid phase does not percolate anymore, large statistical fluctuations occur and the average data depend upon the size of the unit cell because of the classical finite size effect (Fisher, 1971). Finally, the data are shown to verify some upper bounds derived from the best variational techniques. Random packings of monodisperse spheres can be built by sequential deposition of particles (Coelho et al., 1996). The resulting packing was compressed between two plates and the resulting macroscopic coefficients were calculated; the results are compared to the ordered packings. A last interesting result was obtained on a piece of real Fontainebleau sandstone whose pore space was measured by means of Synchrotron computed microtomography (Spanne et af., 1994). As before, the sample is compressed between two plates; the resulting macroscopic coefficients are in good agreement with the ones obtained for reconstructed porous media ; they are found to be larger than the experimental data, which is probably due to the presence of microfissures in the real materials which are not taken into account in the computations as discussed in Poutet et al. (1996). Section 6 is devoted to a discussion of the general properties of the data obtained in this work ; they are summarized in a few graphs as functions of porosity ; they are shown to be smaller than the upper bound of Hashin and Shtrikman (1962) ; an empirical lower bound is found to be the data relative to media derived from site percolation. Some general remarks end up this paper.

2.

POROUS

MEDIA

Let us consider a macroscopic volume V of porous material which is characterized by a length L [cf. Fig. 1(a)] ; it can be separated in a liquid part VL and a solid part V, ; the solid interface is denoted by S,. The external boundary of V is denoted by (! V and various conditions can be applied to it. At some small scale I (the pore scale), the porous medium is assumed to be spatially periodic [cf. Fig. l(b)], i.e. the whole

1590

J. POUTET

et al.

>

X

b

a Fig. 1. Macroscopic

and microscopic

views of a spatially periodic L and 1 are displayed.

porous

medium.

y

=x,-l

The two length scales

medium can be generated by the translation in D independent directions of a unit cell containing pores of arbitrary shapes. D denotes the dimension of space, which may be either 2 or 3. Let I,, I2 and I3 be the basic vectors of the cell; the translations R, leave the porous medium invariant R, = n,I, +nzI, +n,I,,

(1)

where (n,, n2, n3) belongs to Z3. I is a length of the same order of magnitude as the periodic cell. It is assumed that the two length scales 1 and L differ significantly from each other, I<< L or equivalently,

one can define a small parameter ?j = I/L << 1

q (2)

so that the volume V contains a very large number of unit cells. zL, zs, z = zL + ~~ and SP denote the fluid volume, the solid volume, the total volume and the liquid-solid interface inside the unit cell, respectively. The porosity E may be defined as E = SJ(ZL + zs).

(3)

Most unit cells of the porous media analysed in this paper for numerical applications are classical and they have been already described in some of our previous contributions (Lemaitre and Adler, 1990 ; Sallis et al., 1993). These unit cells have the following common features ; they are either square or cubic ; their size is h. Each unit cells is divided into NF elementary squares or cubes of size a. Each elementary domain is filled either with liquid or with solid according to a given recipe. Three families of porous media are considered in this paper. In order to compare our computations to existing results, spatially periodic arrays of spheres were studied first. The second family is composed of deterministic fractals. The first one is the simplest version of the two-dimensional Sierpinski carpets, denoted by SC; starting from a

Random

I

porous

media

1591

Y

-h-

b

C

Fig. 2. (a) Sierpinski carpet : at construction state I (left) and 2 (right). The whole medium be composed of cells identical to these ones. Three-dimensional fractals at construction (b) Fractal foam and (c) Menger sponge. In (b) and (c), the solid phase is denoted

is supposed to stage N = 2. by L,.

partition of the plane by squares of size h, each square is divided into nine subsquares and the centre subsquare is removed from the each initial square. This division and removal process can be continued as many times as desired. The construction stage is denoted by N. Thus the porous medium is spatially periodic at the large scale and fi-xtal inside the unit cell [see Fig. 2(a)]. Two sorts of three-dimensional fractals were studied. The fractal foam (denoted by FF) is a three-dimensional version of the Sierpinski carpet; the space is divided in cubic cells of size h,and each cell is subdivided into 27 elementary cubes ; the central one is then removed. This process can be repeated in order to obtain the medium displayed in Fig. 2(b). The Menger sponge can be constructed in the same way as the fractal foam, but a three-dimensional central cross, made of seven elementary cubes, is removed instead of one central elementary cube only. It should be noticed that both phases of the Menger sponge are connected in contrast with the two previous cases where only one phase was connected. Therefore, two versions of the Menger sponge can be defined ;

J. POUTET

1592

et al.

S-cluster

Z

-N,.a

----+

X Fig. 3. Example

of three dimensional random media derived from site percolation. The boundaries unit cell are indicated by the broken lines. The solid phase is shadowed.

of the

one is composed of interconnected capillaries of various sizes, while the other presents a network of cylindrical obstacles in an unbounded fluid phase. Only the former (MS) is considered in this paper [Fig. 2(c)]. The third family of porous media is made of three classes of random media. The first class is derived from site percolation in two and three dimensions (Fig. 3). Let us consider a cubic lattice where each elementary square or cube of size a is randomly occupied by solid or liquid. The probability of liquid occupation E is the familiar porosity. General properties of these random media have been recalled by Adler (1992) such as the critical probability E, above which the liquid phase is continuous and the correlation length r, which can be roughly defined as the average distance of two sites belonging to the same cluster. Close to percolation, i.e. close to E_ 5 is shown to diverge and is equivalent to

Random porous media
1593

& E E,.

(4)

The values of these quantities were obtained after extensive numerical simulations E, z 0.593,

v = 413. for D = 2,

E, = 0.3117,

v z 0.9,

for D = 3.

(5)

The value of v for D = 2 is an exact result. It is also important to remind that, close to E,, the percolating cluster is fractal and its fractal dimension is about 2.5. Such a random medium derived from site percolation offers an example of a medium which is both fractal and random. The second class consists of random packings of spheres. They are built by simulating the sequential deposition of particles, in a parallelipedic cell with periodicity conditions along the two horizontal directions. After a particle is introduced at a random position above the bed, it settles vertically until it touches a grain already in place. Then, it rolls and glides until it reaches an equilibrium position, where it sits at least on three supporting points. The search algorithm to find this equilibrium position allows all the rotations and all the translations along the bed surface provided that they lower the particle baricenter. The numerical procedure is described in more details by Coelho et al. (1996). Computations were here limited to monodisperse spheres of radius R equal to 4a, where a is the size of the elementary cube ; the size of the unit cell is equal to (2.5~)~; the resulting porosity E is equal to 0.408. Such a packing is illustrated in Fig. 4(a). The third class consists of a piece of real Fontainebleau sandstone which was analysed by means of Synchrotron computed microtomography (Spanne et al., 1994). The major advantage of such a technique is that it provides the local geometry with an excellent resolution ; in the example which is displayed in Fig. 4(b), the size of the voxel is a = 10 pm ; the size of the cell is equal to (30~)~. The porosity of the sample used for the calculations is E = 0.209.

3.

THE EFFECTIVE

ELASTIC PROPERTIES

3.1. General

Consider a porous medium whose solid matrix is perfectly elastic. The equations that govern the elastic behaviour of the material are the basic equations of elastostatics. Let u denote the stress tensor; in the absence of any external volumetric force, the equilibrium equation reads as V*a=O.

(6a)

The strain tensor e is expressed in terms of the displacement d e = (Vd + (Vd)‘)/2,

(6b)

where the superscript t denotes the transposition operator. Only isotropic solid matrices are considered here. Hence, the stress tensor is given by

X

Y

Fig. 4. Examples of random porous media. (a) Random packing of monodisperse spheres obtained by sequential deposition ; the porosity is 0.408. (b) The pore space of real sample of Fontainebleau sandstone obtained with a Synchrotron computed microtomograph; the size of the voxel is (10 pm)‘; the porosity is 0.209. The two samples are squeezed between two horizontal solid plates of thickness equal to a. The lower plate only is displayed.

c;-

Z

Random porous media d =

1595

1tre*I+2pe,

(6~)

where ,I and .Uare the familiar Lame coefficients. These equations have to be supplemented with boundary conditions at the solidpore interface S,. Let us simply assume that no external force is exerted at this interface whose unit normal is n a-n=0

on&.

(64

Sometimes, one needs to apply an external force T at the pore-solid interface; this occurs for instance when the medium is filled by a pressurized liquid. The previous equation is thus replaced by

Equivalently, placement d

the equilibrium

a-n=T

on&.

equation

can be expressed in terms of the dis-

@+p)V(Vsd)+pV2d

(6e)

= 0

(74

together with the condition of no stress at the boundary S, of the solid phase [p(Vd+(Vd)‘)+i(V*d)I]*n

= 0,

c-n = 0

on S,.

(7b)

As mentioned previously, the porous medium is viewed as a spatially periodic porous medium made of an infinity of identical cells. The right theoretical background to derive the macroscopic properties is the theory of homogenization (SanchezPalencia, 1980), which is briefly recalled here. A multiple scale expansion can be applied to the elastostatics equations (6). The major purpose of this expansion is the upscaling in a suitable manner of (6) at some intermediate scale _Y such that 1 <
(8)

Consider any tensorial function f,,(R) that depends upon the parameter II. Due to the existence of the small spatial scale 1 and of the related variations of f,,(R), this function can be regarded as a function of the two variables R and S f,(R) = T(R, S),

(9)

S = R/r].

(IO)

where S denotes the fast variable

Note that generally ‘i must be carefully distinguished from f,, hence the notational change. The spatial derivative may be expressed as Vf, = vJ+

b,T rl

(11)

As suggested by the introduction of the small parameter q, we shall not try to obtain an exact solution to (6), but rather an approximate one within a known error in terms of ‘I. Hence, it is assumed that f, can be expanded in the form

J. POUTET et al.

1596

f, = P(R,S)+rlf’(R,S)+12f2(R,S)

(12)

with the condition that f,(R, S)

periodic in S with the periods (I,, 12, 13).

(13)

By application of (12), the fields CTand d can be expanded as d,(R) = d”(R)+r/d’(R,S)+qZd2(R,S)+

...

a,,(R) = a”(R,S)+qa’(R,S)+q2a2(R,S)+

....

(144 (14b)

The introduction of these expansions into (6) yields a series of equations where each term corresponding to a given order in v]is equated to 0. Let us consider first the expansion in terms of the stress tensor. At the lowest order, it can be easily shown Vs*ao = 0,

n-6’

= 0

on S,.

(15)

The next order yields an equation where terms of both orders are mixed Vs.-a’ +VR-a” = 0,

n-a’

= 0

on S,.

(16)

In order to obtain the macroscopic equation, (16) is integrated over the unit cell ; since a’ is a periodic function of S, one readily obtains VR-(a’)

= 0,

(17)

where the brackets denote the spatial average over the unit cell of total volume z

(18) The local stress a0 can be expressed as a function of the strain. However, because of (1 l), e” is the sum of two terms e” = [VRdo+ (V,d”)‘]/2 + [Vsd’ + (V,d’)‘]/2.

(19)

Hence, with obvious notations a0 = a”(do)+ao(d’).

(20)

This decomposition can be introduced into the local equation (15) ; since do has been assumed to be independent of S [cf. (14a)], (15) reduces to Vs*ao(d’)

= 0

n.a’(d’)

= -n*a”(do)

on S,,.

(21a)

When an external force T has to be taken into account [cf. (6e)], the previous system is modified as Vs*ao(d’)

= 0,

n*a’(d’)

= -n.aO(d’)+T

on S,.

(21b)

The macroscopic stress can be obtained by performing the integration over the unit cell

Random

(a”)

= k

porous

media

(a”(do)+ao(d1))*d3S. sT

1597

(22)

It is time to look at what has been achieved so far. It can be rephrased in the following physical terms according to the terminology of Adler (1992). The general problem consists of two parts. Problem 2 consists of the determination of the effective elastic properties of the porous medium. In order to obtain them, the local displacement d’ should be obtained by solving (21). It is important to notice that (21) is a linear problem ; hence, the solution is a linear function of the “driving force”, which is the unknown displacement do ; more precisely, the solution is a linear function of the strain tensor

E = [V,d” + V,d”)‘]/2.

(23)

It is an easy matter to realize that E is also the macroscopic strain, i.e. the spatial average (18) of the local strains. Thanks to the linearity of (21), one can derive its solution under the general form a’(d’)

= Y:E,

(24)

where 9 is a fourth-order tensor. Subsequent use of (22) and (24) easily yields (a”) = A:E,

(25)

where the fourth-order tensor A represents the effective elastic properties of the porous material. Once A is determined, one can cope with Problem 1, which consists of the determination of do. This can be done by solving (17), together with (25) and the relevant conditions at the boundary aI’ of the macroscopic volume V of the porous medium [see Fig. 1(a)]. This problem will not be further considered here. 3.2. The effective macroscopic properties The material fourth-order tensor A can be simplified whenever the porous medium possesses geometric symmetries ; this is a standard property of solid materials (Landau and Lifshitz, 1970). In general, A has 21 independent components, which can be presented under the abbreviated form

El1 E22 EC.

E33

(26)

E23 E3, El2

where C is a 6 x 6 symmetric matrix whose elements C,, can be numbered in the usual way ; i and j are the line and row numbers, respectively.

J. POUTET

1598

When the porous

medium

is isotropic, (a)

et al.

this relation

can be simplified

into

= 1, trE*I+2p,E,

(27)

where 1, and pL,are the effective Lame coefficients of the medium. There are two additional important cases in view of the following applications to particular media. The first one concerns media with cubic symmetries, i.e. with three mutually perpendicular symmetry planes; an example of such a structure is the Menger sponge displayed in Fig. 2(c). In this situation, one has

c=

Cl’

Cl2 Cl2

0

0

0

c,*

Cl’

Cl2

0

0

0

Cl2

Cl2

Cl,

0

0

0

0

0

0

c&$4

0

0

0

0

0

0

0

0

0

0

c,, 0

(28)

0 c,,

The second one concerns media which possess a translational symmetry along the x-axis ; when the medium is symmetric with respect to the xy- and xz-planes, the tensor C is given by

c=

Cl1

Cl2

Cl2

0

0

0

Cl2

c22

c2,

0

0

0

c12

c2,

c22

0

0

0

0

0

0

c‘$4

0

0

0

0

0

0

0

0

0

0

(29) c,, 0

An example of such a medium is provided by a material yz-plane is a Sierpinski carpet displayed in Fig. 2(a).

0 c,,

whose cross-section

in the

3.3. The numerical routines ; validation and accuracy The problem corresponding to the equation and the boundary conditions (21) is solved via a second-order finite-difference formulation. The mesh spacing A is the same in the three directions and is equal to a fraction a/n of a, the size of the elementary cube. n is at least equal to three, so that the order of the method is equal to two. Equation (21) was discretized by means of the finite volume method in terms of the unknown displacements d’. The external force T may be present or not. Equation (21) was integrated over the elementary volume dlf, which is in most cases a cube of size a/n ; the three unknown components of d’ were taken at the centre of d I’. The divergence theorem was used and quantities such as Vd’ have to be integrated over the faces of d V. When d V is only surrounded by solid volume elements, the approximation of these derivatives is relatively easy, but tedious; three equations with 15 different terms are obtained ; these terms correspond to the unknown displacements at the center of d V and at 14 neighbour elements.

Random porous media

1599

However, the situation is still more complex and one has to take into account the various cases where dV is in contact in one way or another with a liquid neighbour. The amount of work contained in this step is quite large, since 69 different cases are recognized and specifically discretized. Finally, the non-zero elements of the matrix of the corresponding linear system are determined and stored for a given content of the unit cell. The linear system can be symbolized by the equation A*x = b,

(30a)

where x is the vector corresponding to all the unknown displacements. It is solved by using the conjugate gradient method. The explicit writing of the full matrix A is avoided and the scheme proved to be relatively efficient. The convergence test is very simple ; the norm of the rest is estimated at regular iteration intervals ; the program stops whenever the following condition is fulfilled

(job) The precision Ar is of the order of 10P7-10-4. Its influence will be illustrated in Section 4. Typically, the programs were run on an IBM RS6000-560 workstation ; the memory requirement is about 340 (1 -E) (n * NJ3 bytes. The duration of the computations depends on the accuracy to be achieved and on the intricacy of the geometry. For example, highly converged data for arrays of spheres with N, = 20, n = 3 (Z 300,000 unknown displacements) required about 15 h per case. The treatment of random media close to the percolation threshold for the solid phase is particularly time consuming. Once the field d’ is determined, Vd’ is easily derived and the average macroscopic stress (a’) is deduced from (22). For each configuration, one can obtain the tensor A by successively choosing six macroscopic independent macroscopic deformations ; each component of the vector at the right-hand side of (26) is successively set equal to one and the other ones are set to zero. Since the resulting code is complex, it is absolutely necessary to validate it as much as possible before its systematic use. A series of tests was performed that we shall expose now. The first series of tests consists of symmetry tests ; a very efficient configuration in this respect is the three-dimensional spatially periodic cross which corresponds to an infinite cubic array of bars with square cross-section; these bars are parallel to the three coordinate axes and thus mutually perpendicular. This cross is submitted to a constant external pressure; of course, the resulting displacement field is unknown, but it should display symmetries with respect to the three coordinate planes. Another efficient test consists of verifying that the displacements should be left unchanged when the unit cell is shifted with respect to the porous medium. These two tests were successfully passed on the three-dimensional cross. It is known that not many solutions of the elastostic equations exist to which our

computations can be quantitatively compared. Consider an isolated body which is submitted to an external homogeneous and anisotropic force field T

1600

J. POUTET et al.

where n is the external unit normal. The displacement with respect to the centre of gravity of the body is given by

di = -PI+

4P, +P2+Pd

3&2~

1 4

2~’

@lb)

Such a property was verified for a cube and an arbitrary body submitted to a constant pressure field, and finally for a sphere submitted to an artificial anisotropic pressure field. The numerical results were found to be in exact agreement with the analytical predictions. The ratio between the error on the displacement and the maximal displacement is of the order of A,. Systematic numerical errors may be due to three sources, namely the residual A, in the solution of the linear system (30a), the discrete geometry error associated with N, and the discretization error in solving for the displacement field, associated with IZ. The influence of the precision A, on the effective properties was studied for the two particular cases of a Menger sponge at generation 2 and of a random three-dimensional medium of an intermediate porosity E = 0.3, which is close to the percolation threshold E, of the pore space [cf. (5)]. In all these examples, the routine was found to work satisfactorily. The influence of n on the macroscopic coefficients C, was systematically investigated for the two-dimensional Sierpinski carpet at generation 2, which is displayed in Fig. 2(b). n was varied between 3 and 6 for a Poisson’s ratio v = 0.2 and a Young’s modulus E = 1. The computation times required to meet the criterion A, = 10e5 were multiplied by a factor 30. The influence on the coefficients C, was found to be relatively small: all the C, vary by less than 5.10e3 between II = 3 and n = 6; the effective coefficients v2,, and E,, [see (34)] vary by 5.10e3 and 2.5 lop3 respectively. The most precise comparison with former literature results could be achieved for cubic arrays of void spheres, since some attention was devoted to those structures in the past. This opportunity was taken to study the influence of n and of the discretization N,. Cubic lattices of rigid and elastic spheres embedded in an isotropic elastic solid were studied by Nunan and Keller (1984) and by Sangani and Lu (1987) respectively. Nemat-Nasser and Taya (198 1, 1985) and Iwakuma and Nemat-Nasser (1983) studied the influence of periodically distributed voids. Chapman and Higdon (1994) considered cubic arrays of overlapping spheres with an interstitial space filled by a compressible fluid. The influence of n was briefly analysed for a simple cubic close packing of void spheres. For N, = 11 and v = 0.3, n was varied from 3 to 5. C,r and C,2 ranged between 0.713 and 0.715, and 0.230 and 0.231, respectively. Again, the influence of n is limited to the third digit of these coefficients. Our results for cubic arrays of void spheres were compared to the results of Sangani and Lu (1987) and Iwakuma and Nemat-Nasser (1983).Two series of computations with sample sizes N, = 21 and 11 were performed for a Poison ratio v = 0.3 in order

Random

porous

1601

media

Y

-1.4

0.05

I

I

0.1

Fig. 5. Coefficients x(-), porosity E. Data are for:

0.15

0.2

I

I 0.25

0.35

0.3

0.4

0.45

0.5

&

/J(. .), y(---) [cf. (32b)] for cubic arrays of void spheres as functions of the x (Sangani and Lu, 1987), * (Iwakuma and nemat-Nasser, 1983), 0 (our data. N, = 1I), + (our data, N, = 21).

to study the influence of the sphere discretization ; the sphere radius is denoted by R. The discretization rule is the following ; the sphere centre is located at the centre of an elementary cube and any elementary cube whose centre is at a distance less than R is void. The resulting porosity is thus easily computed. The parameters of the computations were chosen as N, = 21,

R = 8.21,

N, = 11,

R = 4.3,

E = E =

The results of these authors were not expressed three coefficients ~1,/?, y defined by [see (25)] A!Jkl =(J+

0.251, 0.256.

(32a)

in terms of the tensor C, but of the

CLY)Giiskl+~L(1+/j)(~rksj~+~ilfi,k)+2Cl(a-B)6,j~,.

Wb)

Our data are compared to these results in Fig. 5. Clearly, they are very close to the results of Sangani and Lu (1987) and Iwakuma and Nemat-Nasser (1983). Moreover, since our precision A, is better than 5.10P4 in both cases, most of the discrepancy is likely due to the sphere discretization, since the agreement is considerably improved when N, changes from 11 to 21. This terminates the preliminary series of checks performed on the routine, which was found to work satisfactorily well. It will be seen in the following sections that the routine satisfies some theoretical properties the two-dimensional and three-dimensional media.

1602

J. POUTET et al.

4.

DETERMINISTIC

STRUCTURES

This section is devoted to the presentation of the results relative to fractal structures. For sake of completeness, it is useful to start by a few properties of materials with symmetry planes.

4.1. Materials

with symmetry

planes

In order to interpret the results, it is convenient to recall the definitions of Young’s modulus and of Poisson’s ratio for isotropic materials. These notions will be extended to materials which possess only symmetry planes. Consider a three-dimensional bar of Section S parallel to the x-axis and submitted to a simple traction F; the bar is made of a homogeneous and isotropic material. The sides of the bar are free to move. The elastostatic equation (6a) can be solved and one obtains the deformation eXX

exx= PIE,

(33)

where p is equal to F/S. The Young’s modulus is given by E = (31+2~)~ I+,Ll

(344

The Poisson’s ratio v is defined as the ratio between the transversal compression eZZ and the longitudinal deformation eXX

(34b) The same calculations can be performed on a 2D solid, i.e. on a ribbon of width I submitted to a force F. The resulting tension is equal top = F/Z. The two-dimensional Young’s modulus and Poisson’s ratio are denoted by EsD and vZD.They are easily calculated as 3L+p 1+2/J

Ezo =4j~c1---==

V

1

E

l-v*’ V

2D=m=jq=

The same traction experiment can be performed on a bar made of anisotropic materials (see Love, 1944). For sake of simplicity, consider first a porous medium with three-mutually perpendicular symmetry planes yz, zx, xy such as the Menger sponge [cf. Fig. 2(c)]. Take a bar made of this material and submit it to a traction F; the cross-section S is a square whose sides are parallel to the y- and z-axes ; again these sides are free to move and the corresponding stresses are equal to zero. It is elementary to show that

porous media

Random

c12

e,, = -

c,,

+

c,

1603

P

e,, =

2 e,,,

c11-

Hence, one can define equivalent

Young’s

Wa)

2c:,

modulus

c,,+c,*

E,, and Poisson’s

ratio vIZ

Wb) This can be extended to materials with only two symmetry planes, but with a translational symmetry along the x-axis such as the Sierpinski carpet [cf. Fig 2(a)]. The tensor C is given by (29). Two traction experiments can be performed. In the first case, the bar axis is parallel to the x-axis. Because of the translational symmetry, it is not difficult to show that the problem is identical to the previous one with the same results given by (35). A second type of experiment can be performed. Take a bar with a square cross-section in the xz-plane ; its axis is parallel to the y-axis ; a traction is exerted along the y-axis and one has to solve for the three deformations eXx,e,, and e,,. The equivalent Young’s modulus Es2 may be expressed as E

c

_

22

-

+

w23-G2)c:2-cI*G

(36a)

22 c,,c22-c;2

Two Poisson’s ratios which is examined exx “2’

=

-<=

can be defined,

CC22

-

c,,c22-c;2

depending

upon

the transversal

e,=

GX12 ’

v23

=

-e,,=c,,c22-cf2.

deformation

c,IG-c:2 (36b)

As before, a two-dimensional Sierpinski carpet can also be considered. Now, a force F is exerted along the y-axis on a ribbon contained in the z-plane. The resulting two-dimensional Young’s modulus E,2,ZD and Poisson ratio v23,2Dmay be expressed as (36~)

(36d) 4.2. Sierpinski carpets The elastostatics equations (21) were solved for the classical ternary Sierpinski carpets displayed in Fig. 2. Since this structure is translationally invariant along the x-axis, computations can be made relatively easily; they were performed for four Poisson’s ratios (v = 0.1,0.2,0.3,0.4) and up to generation 4, whenever possible. For sake of brevity, only the results for v = 0.2 are displayed in Table 1(a). The equivalent Young’s moduli and Poisson’s ratios defined by (35) and (36) were systematically determined. Consider first the case where the material is put under traction along the x-axis.

(4

0.2

V

(cl

0.2

V

@)

0.2

V

0.2

V

0 1 2 3

Generation

0 1 2

Generation

0 1 2 3

Generation

0 I 2 3

Generation

1.1111 0.56 0.3071 0.1776

C,,

1.1111 1.0186 0.9397

C,,

0.5 0.55 0.58

Ratio

0.92 0.92

Ratio

0.99997 0.88889 0.79012 0.70233

C,2

VI2

0.34 0.43 0.49

Ratio

0.88 0.89

Ratio

0.20001 0.2 0.2 0.20001

0.2778 0.094 0.0404 0.0196

C,,

C,, 1.1111 0.79676 0.59157 0.44154

0.2778 0.2432 0.217

0.89 0.89 0.89

Ratio

E,,

Ratio

0.87 0.88 0.88

C,,

1.1111 0.96563 0.8461 0.74355

1. The effective

C,,

0.4167 0.1765 0.079 0.0509

C,,

0.4167 0.3802 0.3472

1 1 1

Ratio

0.72 0.74 0.75

Ratio

0.42 0.45 0.64

Ratio

0.91 0.91

Ratio

1.07406 0.77306 0.57629 0.43171

Ez2

0.4167 0.2606 0.16299 0.10289

C,

1.1112 0.447 0.1984 0.1214

C’“” II

1.1112 1.0036 0.9114

C’fy

0.72 0.75 0.75

Ratio

0.63 0.63 0.63

Ratio

0.40 0.44 0.61

Ratio

0.90 0.91

Ratio

0.20001 0.1661 0.14066 0.11932

V21

0.4167 0.32275 0.25627 0.20455

C,,

0.99997 0.53298 0.29771 0.1737

E,,

0.99997 0.9248 0.8583

E,,

0.83 0.85 0.85

Ratio

0.77 0.79 0.80

Ratio

0.53 0.56 0.58

Ratio

0.92 0.93

Ratio

VIZ

0.82 0.91 0.93

Ratio

0.59 0.67 0.68

Ratio

0.20001 0.14373 0.11626 0.09939

VI2

0.20001 0.19274 0.1876

0.20001 0.16402 0.14953 0.13902

V23

0.2778 0.16255 0.10814 0.07368

C,,

0.72 0.81 0.85

Ratio

0.96 0.97

Ratio

1 0.4815 0.2318 0.1116

EN

1 0.9259 0.857

EN

1.1112 0.80805 0.62068 0.48278

C’“” 22

0.2778 0.19186 0.13994 0.10305

C,,

0.2 0.2 0.2 0.2

vN

0.2 0.2 0.2

vN

0.73 0.77 0.78

Ratio

0.69 0.73 0.74

Ratio

macroscopic properties of the three deterministic fractals. Data are for E = 1, v = 0.2. The ratios of each quantities for successive generations are systematically given. (a) Sierpinski carpets (translational/y invariant along the x-axis) ; (b) fractal foams ; (c) Menger sponges. EN and v,,, in (b) and (c) are the predictions of (49)

Table

a

s 2 z

8

13

g

Random porous media

1605

Obviously, this is a degenerate case where the system undergoes an homogeneous deformation since all the stresses, except o,,, vanish at the boundary of the cell and at the boundary of the inner hole. The coefficients of the carpets can be determined by the following simple argument. If p denotes the tension relative to the surface of the unit cell, the tension which is actually exerted per unit surface of the material is 9p/S, since the solid surface occupies S/9 of the surface of the unit cell. Since the deformation is homogeneous, the Poisson’s ratio of the carpet is equal to the Poisson’s ratio of the material. Hence, for a carpet at generation N, one has E,,

=(I)”

v,, = v.

(37)

This is precisely what is given in Table 1(a) ; this property was also valid for the other values of v. A further comment can be made on this physical situation ; if the porous medium had a translational symmetry along the x-axis and was isotropic around it, one would have [cf. (27) and (28)] c;s; = czj +2&,. This value is systematically compared with CZZin Table l(a). It is seen that the two values are relatively close one to another, as well as the ratios between successive generations. Consider now the second situation where the Sierpinski carpets are put under tension along the y-axis. The two following inequalities are always verified

&z

v,,
forv=0.1,0.2.

vZ3 < v2] < viz

for v = 0.3,0.4.

(39)

The first of these two inequalities means that the deformation is less important when the material is put under tension along the x-axis than along the y-axis; this will be explained simply by means of an approximate renormalization argument. The physical situation seems to be more complex for the Poisson’s ratio. Finally, a last general comment can be made on the power laws. The ratio between two successive generations is systematically given in Table 1(a) for all the computed values which can be expressed as C cz a:,

(40)

where a, is a constant. It is remarkable to see how good (40) is, since it is obtained for most quantities almost immediately within a few percent. However, the constant a, depends upon the quantity under consideration, and upon the Poisson’s ratio. Roughly, a, varies between 0.65 and 1. One may try to derive these results in a simple manner by using real-space-renormalization group transformations ; such transformations were already applied by Bergman and Kantor (1984) for studying the elastic behaviour of regular fractal networks. An elementary renormalization argument consists of viewing the Sierpinski carpet at generation N as made of eight carpets at generation N - 1 which behave almost independently under a tension p = F/S directed along the y-axis. A unit cell is made of two layers with three carpets at generation N- 1 and of one layer with two carpets at generation N- 1. The extreme two layers are submitted approximately to

1606

J. POUTET er al.

the tension p, and the intermediate deformation e; may be estimated as

one to the tension 3p/2. Hence, the overall

Hence, E 22.N

--

6E22,,_,

17

=

Q%E,,,N_ , .

(41b)

Another extreme estimate can be obtained by stating that the middle (N- I)-carpet does not play any role ; in this case, it is obtained E 22.N

=

2E22,N-

I 13

=

Q.667E22.~-

I.

(4lc)

These values can be compared with the values given in Table l(a). The second estimate, which is very crude, gives relatively good results when compared to the numerical data. However, it should be noticed that this type of argument yields a constant Poisson’s coefficient which is very approximate. A possible way to improve the precision of this renormalization argument could be to estimate the ratio between a coefficient at generation N and at generation N- 1 for N = 1, i.e. right at the beginning. It is seen from Table l(a) that this semi-numerical argument is quite successful since such ratios are generally found to be remarkably constant when N increases. Here again, it is found that the asymptotic regime is very quickly reached. Another type of argument can be employed to find these power laws, which has already been successfully used by Thovert et al. (1990) to derive the macroscopic conductivity of regular fractals. When the construction is carried on until N becomes large, the largest void, corresponding to generation 1, sees the surrounding medium as a continuum in a first approximation. This surrounding medium is a Sierpinski carpet at generation N- 1 whose macroscopic properties Cc-’ are supposed to be known. A random medium of properties Ct with a small proportion E of holes of a given shape has macroscopic properties of the form C,,(E) = Cp,(l -k,,,&

(42)

where k,, is a coefficient to be determined. These relations may be used to derive a recursion relation for the Sierpinski carpet whose equivalent porosity is l/9 C;(E) = C;-’ (1 - k,,,/9).

(43)

Such an argument immediately yields a power law. The coefficients k,-,can be derived from a numerical experiment which is separately performed at low porosity E. Results relative to an isolated square hole 2 x 2 located at the centre of a unit cell N, = 20, i.e. at a porosity E = 0.01, are gathered in Table 2. kc,, is readily derived as (44)

2.99

0.992

1.076 1.30 2.01

1.023

1.11 1.35 2.14

0.1 0.2 0.3 0.4

3.16 3.80 6.20

kC’z2

&Z(E)

CP2 CB(E)

0.112 0.267 0.545 0.131

0.114 0.278 0.577 I .43

0.668 0.649 0.578 0.31 I

’

c;,

c%jc;,

4.024 5.517 8.232

1.162

kG 0.871 0.553 0.387 0.085

’

0.455 0.417 0.385 0.357

c:,

0.434 0.040 0.371 0.346

Cd&)

4.41 4.01 3.57 3.12

kc.,.,

carpet at two successive generations

c$j/C&

of the ratio between the ejyective co@cients ofa Sierpinski the dilute approximation (42)

V

Table 2. Determination

0.511 0.554 0.603 0.653

Czlc~‘c’

by means of

z 3 s m

z s

E E

~

1608

J. POUTET et al.

Some of the ratios Cc/Cfl-’ are then deduced from (43) and they can be compared with the “asymptotic” ratios given in Table l(a). The complete numerical results show that the quality of the comparison depends upon the precise value of the Poisson’s ratio. The three coefficients display various behaviours. CFZ/C$‘;’ is better compared to the numerical ratio for small values of v; C&/C!; is larger than the numerical values at small v and smaller at large v ; finally, Cg/Cy; ’ is in agreement with the numerical predictions only for large v. A last check can be performed on the data relative to the Sierpinski carpet by using some recent results valid for composites (cf. Thorpe and Jasiuk, 1992). The Young’s modulus of a two-dimensional material containing holes of a given shape is independent of the Poisson’s ratio of the matrix material. In order to verify this property, one has to use the relation (34~) which relates E2,, to E and v. This property was seen to be extremely well verified. Of course, this property has no predictive value, but can be used as an additional a posteriori check of the precision of the numerical computations. It might be important to note that the renormalization argument (38) also predicts results which are independent of v& ; however, this might be a fortunate coincidence. The average ratio between successive generations is given by ___ EN- I = 0.751.

(45)

2D

The data relative to the Poisson ratio can be rationalized according to the formula valid for two-dimensional composites (see also Day et al., 1992).

(46) Hence, the two-dimensional Poisson’s ratio vfD tends toward a universal value v, irrespective of the value of the matrix Poisson’s ratio. Its variations with N derived from the numerical results were found to be in very good agreement with the predictions of (46) evaluated with the asymptotic value v, = 0.118 ; this comparison holds for v& = 0.1, 0.2, 0.3 and 0.4.

4.3. Fractal foams Fractal foams are presented in Fig. 2(b). Computations were made for two generations only because of the necessary memory space. Moreover, an adequate precision could not be obtained for N = 3. Generally speaking, A* is equal to 10d6 ; for N = 3, only a precision A, = 0.04 could be obtained in a reasonable amount of time. Results were derived for v = 0.1, 0.2, 0.3 and 0.4 ; only the data for v = 0.2 are given in Table l(b). General comments are as follows. Though it is not strictly speaking isotropic, the fractal foam is very close to isotropic in character, since one can define [cf. (27) and

CW

Random

porousmedia

1609

G;) = c,,+2c,,.

(47)

This value is compared to C,, and is seen to be in good agreement with it. Power laws are obtained for the three equivalent coefficients C,,, C,, and C,. Except maybe for C,, the exponent of the power law depends upon the Poisson’s coefficient. Hence, a renormalization argument such as (41) is doomed to fail. A unit cell is made of two layers with nine foams at generation N- 1 and one intermediate layer with eight such foams. The extreme two layers are approximately submitted to the tension p, and the intermediate one to the tension 9p/8. The overall deformation ez can be derived as e:=3 1(

1 EN_

9 + 8E,,

1

(484

+ E2?,Np, 1 p.

Hence, EN = 24E,_ ,125 = 0.96E,_, .

(48b)

When the intermediate layer is assumed not to play any role, one readily obtains EN = 8E,,_, 19 = 0.89E,_,

(48~)

These figures can be compared to the actual values given in Table l(b). Only the order of magnitude is correct. Let us now try the renormalization argument based on an expansion at low porosity similar to (42). The actual value of the porosity in (42) is l/27 which is quite low. Such an argument may be successful since the two values of the ratios between the successive generations are remarkably close. One can be more specific and use the expansion which was obtained by Dewey (1947) for a dilute suspension of spherical inclusions. Hence, when a fractal foam at generation Nis viewed as a dilute suspension of porosity E = l/27 of spherical voids in a medium of macroscopic properties corresponding to generation N - 1, the following relations hold for the Young’s modulus and the Poisson’s ratio 3 (l-v,~,)(9+5v,+,) l-2-7_5v,m, 3 (1 -vi-,)(1

It

vN

=

vN--l

-5v,-,) E

1+~(7-5v,,lvN_,

1, 1

c

.

should be noticed that whatever the initial value of vo, the Poisson’s ratio tends toward vr‘ = 0.2. Hence, the renormalization argument predicts a universal asymptotic regime for any value of the Poisson’s ratio EN = 25EN_ ,127,

vN -+ 0.2,

N -+ 0~.

(50)

The complete results show that the agreement between the numerical and the analytical value of the Young modulus is excellent. However, it should be noticed that the renormalization model predicts the following trends for the Poisson’s ratio which depends upon the initial value of v

1610

J. POUTET et al. vg < 0.2, vg

=

0.2,

vg > 0.2,

vN

t 0.2,

v&I= 0.2, VN 10.2.

(51)

Moreover, v,~ is always found to be a slowly decreasing function of the generation N. Hence, the agreement for vNis less good than for EN. 4.4. Menger sponges

The last classical example of deterministic fractal is the Menger sponge displayed in Fig. 2(c) ; its properties were also systematically determined for v = 0.1, 0.2, 0.3, 0.4. The same set of comments as for the fractal foam could be made about the results for v = 0.2 which are gathered in Table l(c). The isotropic character of the sponge can be checked by calculating Cyl, according to (47) ; this approximation is obviously worse than before. Power laws are certainly obtained in the limit N+ co, but the ratios between successive generations are not as constant as they were for fractal foams. However, it is remarkable that the Young’s modulus is only slightly dependent upon the Poisson’s ratio. This is encouraging for a simple renormalization argument similar to (41). A unit cell is made of two layers with eight sponges at generation N- 1 and one layer with four sponges. The extreme two layers are approximately submitted to the tension 9p/8 and the intermediate one to the tension 9p/4. Hence, the overall deformation eZZ can be derived as in (41). The asymptotic behaviour of the Young’s modulus is EN = 2EN_,/3.

(52)

This estimation is not far from the numerical findings. The dilute approximation (49) was tentatively tried. Of course, the inclusions are not spherical anymore and the porosity is equal to 7/27 which is not small when compared to 1. It was found that the predictions of the Young’s modulus are better for small values of the Poisson’s ratio v. The comparison for the Poisson’s ratio vNis not good ; as for the fractal foams, v,~ is always found to be a decreasing function of N, irrespective of the initial value v [cf. (50)].

5.

RANDOM MEDIA

Three types of random media were analysed, namely random media derived from site percolation, random packings of spheres and a sample of real Fontainebleau sandstone. 5.1. Site percolation

Examples of three-dimensional random structures are displayed in Fig. 3. Most of the computed results are gathered in Figs 6 and 7. Detailed results relative to the average porosity E = 0.3 were analyzed. For the seven samples which were generated with N, = 10 the actual porosity fluctuates by

Random

porous

1611

media

05.

0.4 -

0.1

"0

0.2

0.4

0.3

V”

0.5

(b)

-

0.2 0.16 -

5+

1

a-

0.08 0.06

01 0

1 0.1

0.2

0.3

0.4

0.5

0.6

&

Fig, 6. Three-dimensional random media derived from site percolation. The macroscopic Young’s modulus /? and the macroscopic Poisson’s ratio v* for E = 1 and v = 0.2 are given as functions of the porosity E in (a) and (b), respectively. The horizontal and vertical bars indicate the statistical fluctuations of the porosity and of the macroscopic coefficients. The size NC of the unit cell is given by a number (A’, = 5, 8, 10 or 12). The solid line corresponds to the estimation (56).

J. POUTET

1612

(a)

06

et al.

k*

0.5 0.4 0.3 0.2 0.1 0 0

01

0.2

0.3

0.4

0.5

0.6

0.7

0.1

Q2

03

04

0.5

0.6

0.7

Ye b)

0.4 0.3 0.2 0.1 0 0

Fig. 7. (a) The average bulk modulus K* and (b) shear modulus p* as functions of the porosity E for random media. Data are for: E = 1, v = 0.2, NC = 5 (a), 8 (A), 10 (0). The solid lines are the HashinShtrikman upper bound (57) ; the broken line is Miller’s upper bound (58) with G = 0.164.

+ 7% and the effective mechanical properties by + 14%, which is twice the previous value. The samples of porous media which are generated by site percolation are isotropic only in an average sense. Each sample is anisotropic ; however, it is interesting to compare its properties to the ones of a perfectly isotropic sample, for which one should have because of (27) c,, = c13, c,, = c,, +2c44 = Cl3 f2C44.

(53)

These relations were found to be well verified. Systematic calculations of the effective properties as functions of the porosity and of the size NC of the unit cell were performed ; NC is ranging from 5 to 12 ; the number of samples generated for each porosity is of the order of 10. In order to save computational time only one macroscopic strain tensor

Random porous media

1613

was used. The coefficients C,, and C,, are equal for an isotropic medium ; their average value CyV,can be used to calculate the equivalent Young modulus E* and Poisson ratio v* [cf. (27) and (35)] (54a)

(54b) (54c) These quantities were computed for each configuration and then averaged. It turns out that if these two operations are performed in the reverse order, the final results are identical within three digits. The numerical results relative to E* and v* are displayed in Fig. 6. It is interesting to focus our attention to the region which corresponds to the critical porosity E, = l-e, = 0.688 [cf. (5)] where the solid matrix stops percolating. First, the statistical fluctuations should become very large ; this is more clearly visible on v* which does not tend toward 0 when E + EL.Second, the results depend upon the size of the unit cell according to the classical finite-size effect (Fisher, 1971) ; numerically, E* is found to be a decreasing function of N,, and v* an increasing one (cf. Fig. 6). Close to EL,E* is expected to vary as E* = k(~;--E)~.

(55)

The critical exponent z should be estimated by the finite-size technique (cf. Fisher, 1971), i.e. by computing E* as a function of N, at E = e:. Because the computations become too long, this could not be done. Instead, the following estimation was based on all the results, by a least square fit with the constraints E*(O) = 1 and E*(E,) = 0 (56) This result compares successfully with the numerical data in Fig. 6(a) for E < EL.It is difficult to avoid noticing that r is close to the value of the conductivity exponent whose value is currently estimated as t = 1.92 (Clerc et al., 1990). However, the overall fit (56) is not meant to be a determination of z, which would require a much more careful and expensive finite-size effects analysis. Because of this imprecision, the exponent 1.74 is outside the widely accepted bounds 1 +Dv < T < +2v [Kantor and Webman (1985), Sahimi (1986), Roux (1986)] ; these variational bounds were indeed verified by the numerical simulations on discrete bond percolation networks by Zabolitzky et al. (1986) and by Limat (1988).

1614

J. POUTET et al.

The raise of v* close to E: [Fig. 6(b)] is consistent with the conjecture of Bergman and Kantor (1984) that K*/,u* = D/4 near the threshold (i.e. v* = 0.2 for D = 3). Finally, it was found useful to compare our results to the bounds derived by the variational approach. This comparison is made in Fig. 7 in terms of the bulk modulus K* (which is equal to A*+:@*) and of the shear modulus p*. Since the voids have zero mechanical properties, the lower bounds are trivially zero. The first bounds are due to Hashin and Shtrikman (1963) as further generalized by Walpole (1966). When only Eis known, K* is the best possible bound for all porosities ; p* is the best bound only in the limit F + 0. K* and p* verify (57a)

Another useful bound is due to Miller (1969). Since the space is isotropically and homogeneously paved by elementary solid or liquid cubes in a totally uncorrelated manner, these random media belong to the general class of the “symmetric cell materials” ; note that the term “cell” refers here to the elementary cubes used to pave the space, rather than to the constitutive pattern of the unit cell. The best upper bound when it only involves the porosity E and the shape of the cells is l--E 3K E ‘+ z 3-4~+9(2e-

(58) 1)G)‘

where G has been estimated as 0.164 (cf. Thovert et al., 1990) Our data are found to always satisfy these upper bounds. 5.2. Random packings of spheres Some additional explanations are needed to describe the computations that were performed on the random packing of spheres. The packing is generated according to a continuous scheme, but it is subsequently discretized in elementary cubes of size a for the mechanical calculations; the radius R of the spheres is 4~2,and the size of the unit cell is (25~)~. The packing is isotropic in the xy-plane, because of its construction. Two solid plates of thickness a are added on top and below the packing, as displayed in Fig. 4(a) and a uniform pressure equal to 1 is externally applied on them. Remember that the packing is spatially periodic in the xy-plane. The material properties were chosen to be E = 1 and v = 0.2. Computations were extremely slow to converge, though the requirement on the precision was lowered down to A, = 0.045. This is easily understandable, since the contacts between the spheres are limited to a cube; thus, it is difficult with a uniform

Random porous media

mesh as used here to obtain an outstanding of the stress were calculated as fsI = &LY>

1615

precision.

The two average components

(59a)

+ Co.,,.)),

(59b)

Gil = (a,,>. The deformation

along the z-axis is calculated

as

e,i = (e,,).

(59c)

The porous medium resulting from deposition is assumed to be isotropic; this simple hypothesis enables us to calculate the macroscopic coefficients with the help of (27) which reduces to CL = J.,e,,,

011= (n, +2/L,&

(6’3)

A, and 11,can be deduced from these two relations. Finally, the macroscopic coefficients E* and v* are derived from I, and pL,by (34) which is valid for isotropic materials. When these relations are applied, the following values are obtained I,, = 0.0585,

pL,= 0.123,

E* = 0.286,

v* = 0.161,

These values can be compared to the ones obtained they are not found to be drastically different.

K* = 0.141.

(61)

for cubic arrays of spheres and

5.3. Geological media This example is of particular interest because it is based on the real piece of material displayed in Fig. 4(b). The size of the cell is equal to (30~)‘. where a = 10 pm is recalled to be the size of the model. The computations are made along the same methodology as for random packings of spheres, i.e. the medium is squeezed between two horizontal plates on which a unit force is exerted. The macroscopic properties are deduced with the help of (59) and (60). They are given by i, = 0.149,

pC = 0.214,

E* = 0.645,

Similar sets of results are available from responding values for a porosity of 0.207 are (i) reconstructed

Poutet

K* = 0.332.

et ul. (1996)

(62a)

where

the cor-

K* = 0.27.

(62b)

media :

2, = 0.127, (ii) experimental

V* = 0.167,

/i, = 0.228,

E* = 0.54,

I!* = 0.179,

: E* = 0.45,

K* = 0.15.

(62~)

The comparison of these three sets of numbers is interesting. The numerical data obtained either for the reconstructed media or for the real sample are close one to another; however, they are significantly larger than the experimental data. The discrepancy with the experiments seems to be due to the presence of microcracks which weaken the real material as discussed by Poutet et al. (1996).

1616

J. POUTET

et al.

The Young modulus calculated for the real sample is larger for the reconstructed material ; however, one should be careful general conclusions should be drawn from a single datum. If it again the reconstruction process yields a good approximation this general statement will be corrected below by comparison derived from site percolation.

than the one derived at this point since no can be said that once of the real material, with random media

0.10.08-

x

0.060.040.02.

E 0 0 0.1 0.2 0.3 0.7 0.4 0.5 0.6 Fig. 8. (a) The effective Young’s modulus E*/E, (b) the effective Poisson ratio Y*, (c) the effective bulk modulus K*/& and (d) the effective shear modulus ~~,ih as functions of the porosity E for all media. I?,, = 1; Y,,= 0.2, The line (.-.-.) is the Hashin-Shtrikman’s upper bound. Data are for: site percolation random media (O), reconstructed Fontainebleau sandstones (O), real Fontainebleau sandstone (@), three-dimensional fractals (*), full or hollow spheres in a cubic close packing ( x ), random packing of spheres (8).

1617

Random porous media

6.

DISCUSSION

AND CONCLUDING

REMARKS

All the results obtained for the previous structures for E = 1 and v = 0.2 are gathered in Fig. 8 where they are displayed as functions of porosity. The data from Poutet et al. (1996) for reconstructed media, enriched by three new samples for each porosity, are also included. The following general comments can be made. p* and K* are compared to the upper Hashin and Shtrikman’s bound (1962) and they are found to lie below this limit ; it should be noticed that the largest values are obtained for the real Fontainebleau sandstones and for the fractal three-dimensional

0.7 0.6 0.5 0.4 0.3 0.2 0.1 O0

Fig. 8.-Continued.

1618

J. POUTET

et al.

structures. Moreover, as a matter of fact, an empirical lower bound to most results is provided by the data relative to the random media based on site percolation. Another interesting feature is that the correlation function does not play an important role on the mechanical properties in the range of parameters that is studied here ; the data for the site percolation and the reconstructed media are indeed very close. This feature was not observed before on other types of transport where the influence of the correlation is very important; for instance, permeability was addressed by Adler et al. (1990) and by Lemaitre and Adler (1990) ; a possible explanation is that the solid concentration is rather large and that the detailed structure does not matter in this range. It was not found useful to display on this figure the experimental data given in Poutet et al. (1996), since they are significantly lower than the numerical ones, most probably because of the presence of microcracks. Results for the sphere packings do not differ very much one from another; this might be due to the fact that the porosities are very similar for ordered and disordered packings, and for full and empty spheres. Finally, it is remarkable to note that the results relative to three-dimensional fractals (i.e. Menger sponges and fractal foams) appears to lie on a single curve, when plotted as functions of concentrations. These general remarks terminate this paper. One can conclude by saying that the determination of the macroscopic mechanical properties of most structures of fundamental or applied interest is now possible, thanks to the recent developments of the computational possibilities. We have addressed the properties of several porous media whose discrete versions have stimulated a lot of interest since they provided novel ways to describe the geometry of real materials. Like many of their properties, the elastic coefficients of fractal media are ruled by power laws. These can often be derived, at least approximately, by use of simple renormalization arguments. The Young modulus of three-dimensional random media follows also a power law, over the whole range of porosity. Whenever possible, these numerical results were compared to experimental data.

ACKNOWLEDGEMENTS Most computations were performed at CNUSC is gratefully acknowledged.

(subsidized

by the MESR)

whose

support

REFERENCES Adler, P. M. (1992) Porous media : Geometry and Transports. Butterworth/Heinemann. Adler, P. M., Jacquin, C. G. and Quiblier, J. A. (1990) Flow in simulated porous media. ht. J. Multiphase Flow 16, 691-712. Bergman, D. J. and Kantor, Y. (1984) Critical properties of an elastic fractal. Phys. Rev. Lett. 53,511-514. Bout&a, M. J. (1992) Elements of poroelasticity for reservoir engineering. Rev. Inst. Fran. P&role 47,479-490. Bout&a, M. J., Piau, J.-M., Kessler, N., Boisson, M., Bary, D. and Fourmaintraux, D. (1994) SPE/ISRM 28093, SPE/ISRM Symposium on Rock Mechanics, Delt 1994.

Random porous media

1619

Chapman, A. M. and Higdon, J. J. L. (1994) Effective elastic properties for a periodic bicontinuous porous medium. J. Mech. Phys. Solids 42,283-305. Cleary, M. P., Chen, I.-W. and Lee, S.-M. (1980) Self-consistent techniques for heterogeneous media. J. Engng Mech., ASCE 106,861-887. Clerc, J. P., Giraud, G., Laugier, J. M. and Luck, J. M. (1990) The electrical conductivity of binary disordered systems, percolation clusters, fractals and related models. Adz,. Phys 39, 191-309. Coelho, D., Thovert, J.-F. and Adler, P. M. (1996) Geometrical and transport properties of random packings of spheres and anisotropic particles, submitted. Coussy, 0. (1991) Mecanique des milieux poreux, Ed. Technip. Day, A. R., Snyder, K. A., Garbo& E. J. and Thorpe, M. F. (1992) The elastic moduli of a sheet containing circular holes. J. Mech. Phys. Solids 40, 103 l-l 05 1. Dewey, J. M. (1947) The elastic constants of materials loaded with non-rigid filler. J. A&. Phys. 18, 578. Fisher, M. E. (197 1) The theory of critical singularities. Critical Phenomena, Proc. 5 I st Fermi School, Varenna, Italy (ed. M. S. Green) Academic Press, New York. Garboczi, E. J. and Day, A. R. (1995) An algorithm for computing the effective linear elastic properties of heterogeneous materials : 3-D results for composites with equal phase Poisson ratios. J. Mech. Phvs. Solids 43, 1349-1362. Hashin. Z. and Shtrikman, S. (1962a) On some variational principles in anisotropic and non homogeneous 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 elastic behavior of multiphase materials. J. Mech. Phys. Solids 11, 127-140. lwakuma, T. and Nemat-Nasser, S. (1983) Composites with periodic microstructure. Compuws Sirucf. 16, 13-20. Joshi, M. (1974) Ph.D. Thesis, University of Kansas. Kantor, Y. and Webman, I. (1984) Elastic properties of random percolating systems. Phw. Rev. Lett. 52, 18914. Laurent, J., Bout&a, M. HJ., Sarda, J.-P. and Bary, D. (1993) Pore-pressure influence in the poroelastic behavior of rocks : experimental studies and results. SPE Formarion Evaluation 117-122. Lemaitre, R. and Adler, P. M. (1990) Fractal porous media. IV-Three-dimensional Stokes flow through random media and regular fractals. Trunsp. Porous Mediu 5, 3255340. Limat, L. (1988) Rotationally invariant elasticity in a planar network. Phys. Rer. B. 38, 512519. Love, A. E. H. (1944) A Treatise on the Mathematical Theory of Elusticitv. Dover, New York. Miller, M. N. (1969) Bounds for effective bulk modulus of heterogeneous materials. J. Math. Phys. lo,200552013. Nemat-Nasser. S. and Taya M. (1981) Body containing periodically distributed voids. Q. Appl. Matlz. 39,43-59. Nemat-Nasser, S. and Taya, M. (1985) On effective elastic moduli of an elastic body containing periodically distributed voids, comments and corrections. Q. Appl. Ma& 43, 187-188. Nunan, K. C. and Keller, J. B. (1984) Effective elasticity tensor of a periodic composite. J. Mech. Phys. Solids 32, 259-280. Phan-Thien, N. Effective lame’s moduli of cubic arrays of elastic spheres embedded in an elastic matrix. J. Appl. Math. Phys. (ZAMP) 34, 387-397. Poutet. J., Manzoni, D., Hage-Chehade, F., Jacquin, C. J., Bout&a, M. J., Thovert. J.-F. and Adler, P. M. (1996) The effective mechanical properties of reconstructed porous media. Int. J. Rock Mech. Min. Sri. 33,409415. Quiblier, J. A. (1984) A new three-dimensional modeling technique for studying porous media. J. Call. Interj: Sci. 98, 84-102. Roux, S. (1986) Relation between elastic and scalar transport exponents in percolation. J. Phys. A. 19, L351-L356.

1620

J. POUTET et al.

Russel, W. B. and Acrivos A. (1972) On the effective moduli of composite materials : slender rigid inclusions at dilute concentrations. ZAMP 23,434464. Sahimi, M. (1986) Relation between the critical exponent of elastic percolation networks and the conductivity and geometrical exponents. J. Phys. C. 19, L79-L83. Sallb, J., Thovert, J.-F., Delannay, R., Prevors, L., Auriault, J.-L. and Adler, P. M. (1993) Taylor dispersion in porous media. Determination of the dispersion tensor. Phys. Fluids A5, 2348-2376. Sanchez-Palencia, E. (1980) Non-Homogeneous Media and Vibration Theory. Springer-Verlag, Berlin. 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. Snyder, K. A., Garboczi, E. J. and Day, A. R. (1992) The elastic moduli of random twodimensional composites : computer simulation and effective medium theory. J. Appl. Phys. 72,5948-5955. Spanne, P., Thovert, J.-F., Jacquin, C. G., Lindquist, W. B., Jones, K. W. and Adler, P. M. (1994) Synchrotron computed microtomography of porous media : topology and transports. Phys. Rev. Lett. 73, 2001-2004. Stauffer, D. (1985) Introduction to Percolation Theory. Taylor and Francis, London. Thorpe, M. F. and Jasiuk, I. (1992) New results in the theory of elasticity for two-dimensional composites. Proc. R. Sot. Lond. A438, 531-544. Thovert J.-F., Wary, F. and Adler, P. M. (1990) Thermal conductivity of random media and regular fractals. J. Appl. Phys. 68, 3872-3883. Vincke, 0. and Bout&a, M. J. (1992) An estimation of terminal bulk moduli of sandstones using their petrographic and petrophysics description. Structure et comportement mecanique des geomateriaux, Colloque RenC Houpert, September 1992, Nancy. Walpole, L. J. (1966) On bounds for the overall elastic modulus of inhomogeneous systems I. J. Mech. Phys. Solids 14, 151-162. Watt, J. P., Davies, G. F. and O’Connell, R. J. (1976) The elastic properties of composite materials. Rev. Geophys. Space Phys. 14, 541-563. Yao, J., Frykman, P., Kalaydjian, F., Thovert, J.-F. and Adler, P. M. (1993) High-order moments of the phase function for real and reconstructed model porous media : a comparison. J. Colloid Interface Sci. 156,478490. Zabolitzky, J. G., Bergman, D. J. and Stauffer, D. (1986) Precision calculation of elasticity per percolation. J. Stat. Phys. 44,21 l-223.