## Abstract

Most studies of effective properties of random heterogeneous materials are based on the assumption of the existence of a representative volume element (RVE), without quantitatively specifying its size *L* relative to that of the micro-heterogeneity *d*. In this paper, we study the finite-size scaling trend to RVE of the Darcy law for Stokesian flow in random porous media, without invoking any periodic structure assumptions, but only assuming the microstructure's statistics to be spatially homogeneous and ergodic. By analogy to the existing methodology in thermomechanics of random materials, we first formulate a Hill–Mandel condition for the Darcy flow velocity and pressure gradient fields. This dictates uniform Neumann and Dirichlet boundary conditions, which, with the help of two variational principles, lead to scale-dependent hierarchies on effective (RVE level) permeability. To quantitatively assess the scaling trend towards the RVE, these hierarchies are computed for various porosities of random disc systems, where the disc centres are generated by a planar hard-core Poisson point field. Overall, it turns out that the higher is the density of random discs—or, equivalently, the narrower are the micro-channels in the system—the smaller is the size of RVE pertaining to the Darcy law.

## 1. Introduction

Most studies of random heterogeneous materials are based on the assumption of the existence of a representative volume element (RVE), and are focused on the properties of the RVE without quantitatively specifying its size *L* relative to that of the micro-heterogeneity *d* (e.g. Berryman & Milton 1985; Nemat-Nasser & Hori 1999; Milton 2002; Torquato 2002). More specifically, the RVE is taken to satisfy the so-called ‘separation of scales’,(1.1)where *d* is the *microscale* (average size of heterogeneity, figure 1*a*), *L*_{Macro} is the *macroscale* of the continuum problem considered. Thus, the right-hand side inequality in (1.1) ensures the applicability of continuum mechanics with Darcy's law as the constitutive relation governing the response of the RVE. However, the left-hand side inequality in (1.1) is postulated rather vaguely (and sometimes involving ), without clearly specifying what is meant by , or, sometimes, postulating a periodic unit cell. In this paper, we study this inequality for flow in porous media within a general framework of scale-dependent homogenization of random media (Ostoja-Starzewski 2001) without any periodicity assumption, in the vein of Hill's (1963) definition of RVE.

In order to attain the said homogenization, two conditions must be satisfied: (i) the microstructure's statistics must be spatially homogeneous and ergodic; (ii) the material's effective constitutive response must be the same under uniform boundary conditions of Dirichlet (essential) and Neumann (natural) types. This approach has been applied in studies of various effective continuum properties in solid mechanics—examples are: classical elasticity (e.g. Huet 1990; Sab 1992; Kanit *et al*. 2003), thermal conductivity (Ostoja-Starzewski & Schulte 1996), viscoelasticity (Huet 1995), physically nonlinear elasticity (Hazanov 1999), elasto-plasticity and rigid-plasticity (Jiang *et al*. 2001; Ostoja-Starzewski 2005; Li & Ostoja-Starzewski in press), linear thermoelasticity (Du & Ostoja-Starzewski in press) and finite (thermo)elasticity (Khisaeva & Ostoja-Starzewski 2006).

In this paper, we consider a steady, single-phase Stokesian fluid flow in random porous media with explicit account of a spatially disordered microstructure. In phenomenological, deterministic continuum mechanics, such a flow is described by Darcy's law (1856). With the slow flow in a porous medium, also considered as incompressible and viscous, Darcy's law is governed by the equations (Dullien 1979)(1.2)(1.3)where is the Darcy (volume averaged) velocity, is the applied pressure gradient driving the flow, is the fluid viscosity and is the permeability, a second-rank tensor that depends on the microstructure of the porous medium. Equation (1.3) is the continuity equation, wherein *g*(*x*) is the source/sink term. The tensor is usually considered a material property, since it definitely depends on the porous microstructure. (In this paper, we interchangeably use the symbolic notation for first- and second-rank tensors and or indicial notations and , whichever is more convenient.)

Clearly, equations (1.2) and (1.3) show that Darcy's law shares the same mathematical formulation as conductivity problems. However, as stated by Torquato (2002, p. 348), ‘As in the trapping problem, there is no local constitutive relation in flow problem. Furthermore, both fluid permeability and trapping constant are scale-dependent properties, in contrast to both the effective conductivity and elastic moduli, which are scale-invariant properties.’ As is well known, the local governing equation of fluid in pore spaces is the Stokes equation.

At this point, let us briefly review prior work on permeability problems in porous media. Evidently, due to the complex, randomly structured nature of such media, deterministic solutions are not possible. Statistical methods are therefore often used for the determination of rigorous (upper/lower) bounds on the effective (RVE level) permeability. In the very early work done by Prager (1961), a variational principle of energy dissipation rate was applied and, with a proposed trial stress field, the upper bound was obtained. More rigorous bounds have been proposed by Berryman & Milton (1985) and a different trial stress with two point bounds was used. They were able to reduce the upper bound by a factor from 9/10 to 2/3. Another important approach has been invented and developed in Rubinstein & Torquato (1989). In their derivation of Darcy's Law from the Stokes equation, an energy-formed functional of permeability was formulated with the ensemble of a ‘tensor velocity field’ and a ‘vector pressure field’. The variational principle was applied on this energy-formed functional, and both upper and lower bounds were obtained with properly selected trail functions. While all the aforementioned works gave bounds on effective permeability, the size of RVE was only required to be much larger than the microscale.

In our work, we introduce a mesoscale domain of size *L*, recall equation (1.1), and use a non-dimensional parameter(1.4)The response of that domain is therefore a function of (i) mesoscale *δ*, (ii) the boundary conditions and (iii) the actual microstructure, such as that shown in figure 1*a*. Since the microstructure, in the ensemble sense, is random, the mesoscale response is statistical, but tends to become deterministic as *δ* grows. It is this tendency in function of *δ* which is the focus of our study.

## 2. Mathematical formulation

### (a) Variational principles

In the general studies on property bounds, the variational principles usually provide the basic starting point (Prager 1961; Huet 1990; Jeulin & Ostoja-Starzewski 2001; Torquato 2002). Now, as mentioned before, Darcy's law (1.2) and (1.3) is not a local constitutive relation, since the slow flow in voids is governed by the Stokes equation:(2.1)where is the actual fluid velocity. In any chosen element d*V* of the porous medium, the energy dissipation rate *ϕ* in a fluid of viscosity *μ* is usually calculated as(2.2)where is void region, is a symmetric, zero-trace stress tensor, which can be obtained by solving the Stokes equation directly.

Let us now consider a trial stress field, , satisfying(2.3)where is a scalar function, a trial pressure field, and(2.4)where is a prescribed pressure on , the surface of the domain *B* of the porous medium. Prager (1961) and Berryman & Milton (1985) stated the variational principle for an element as follows:(2.5)where is the trial dissipation rate calculated with the trial stress function. As we shall see in §2*c*, this gives the upper bound on the effective permeability according to Prager (1961). The equality occurs only when is equal to the true stress, and this is equivalent to the law of energy conservation. This fact will be used in the subsequent numerical simulations.

On the other hand, the variational principle based on Darcy's law (1.2) and (1.3) can be expressed in terms of the pressure gradient and the Darcy velocity. Hence, the energy dissipation rate in the porous domain can also be written as(2.6)Here, we use another form of Darcy's law, in which the pressure gradient is a function of Darcy's velocity:(2.7)with *R*_{ij} being the resistance tensor of the porous medium.

For the effective tensors—that is, on the scale *L* where the separation of scales (1.1) holds—we have the following relation:(2.8)where I is the second-rank identity tensor. However, when the element scale is smaller than the RVE size, the foregoing relation is generally not true.

According to the preceding discussion, the energy dissipation can be expressed in two different quadratic forms:(2.9)

(2.10)

If we look at Darcy's law (1.2) and (1.3), the following equation can easily be derived in the absence of sources and sinks inside the material domain:(2.11)Of course, when considering homogeneous and isotropic media, the foregoing equation becomes the Laplace equation.

Clearly, we see an analogy of the ansatz of flow through random porous media to that of mechanics of random solids (Jeulin & Ostoja-Starzewski 2001). That is, the pair (pressure field, Darcy velocity) corresponds to the pair (strain, stress) in anti-plane elasticity, or the pair (temperature, heat flux) in heat conductivity. Then, is treated as the potential energy, while is the complementary energy.

Now, recall the following two variational principles in terms of potential and complementary energies.

Among the admissible solutions for the pressure field in a porous medium domain, the true solution possesses the minimum potential energy.

Among the admissible solutions for the Darcy velocity field in a porous medium domain, the true solution possesses the minimum complementary energy.

As noted earlier, there is another energy formulation with associated variational principles (Torquato 2002). However, as they do not involve boundary conditions, they will not be used in the present study.

### (b) Hill condition in porous media

When considering a mesoscale domain of the heterogeneous medium, the distribution of surface value will be used to evaluate the constitutive response. In the context of elasticity as defined by Huet (1990, 1995), this results in the so-called apparent properties. As a first step, the Hill condition must be satisfied. In an elastic heterogeneous material, that condition can be expressed as follows:(2.12)where is the unit normal vector of a surface element d*S*, is the coordinate of d*S*. Combining this with the average strain and stress theorems in elasticity, this ensures that the energetic definition is equivalent to the mechanical definition of Hooke's law. Hence, the variational principles based on energy can be applied to evaluate the constitutive response of heterogeneous material in terms of the boundary values.

Analogous to this, a similar theory must also be established in permeability of porous media. Thus, in a given porous medium domain *V*, the pressure gradient and the Darcy velocity can be expressed in the form of average values and fluctuations,(2.13)(2.14)The energy dissipation is calculated as(2.15)From the foregoing equation it appears that, in general, the energy dissipation cannot be calculated from surface values, except only when . By inspiration of Hill's condition in elasticity (Hill 1963), the following derivation applies:(2.16)Owing to the continuity equation , we get(2.17)Applying the Green–Gauss theorem, we arrive at(2.18)Therefore, in order to satisfy , the term (2.18) must vanish. We conclude that the Hill condition for Darcy's law is(2.19)To satisfy (2.19), there are three possible boundary conditions:(2.20)(2.21)(2.22)Following the terminology of thermomechanics of solid random media, (2.20) is called a *uniform Dirichlet* (or *essential*) boundary condition, (2.21) is called a *uniform Neumann* (or *natural*) boundary condition, while (2.22) is called a *uniform orthogonal-mixed* boundary condition. The Hill condition (2.19) ensures that the energy dissipation calculated from the surface values of the mesoscale element under either of these conditions is equal to the energy dissipation obtained from the volume integration over that element.

At this point, we introduce apparent dissipation functionals on the mesoscale and , which are homogeneous functions of degree 2 in the volume average pressure gradient and the volume average Darcy velocity , respectively. Thus, we have (Ziegler 1983)(2.23)With this we have a full analogy between the Stokesian flow in porous media and the anti-plane elasticity in a matrix of finite stiffness, containing infinitely stiff inclusions.

Note that there is no body force considered in the above derivations. Also, the above derivations use the tensor index notation for clarity, while further derivations will mostly use the symbolic notation.

### (c) Scale-dependent hierarchies of bounds

We take the porous medium as a set of realizations , defined over the sample space . Thus, for an elementary event , we have a realization , theoretically of infinite extent, from which we cut out a mesoscale specimen , of finite size (figure 1*a*). As we shall see below, that specimen is characterized by apparent (or mesoscale) permeability, , and resistance, , tensors, corresponding to (2.20) and (2.21), respectively, according to(2.24)

(2.25)

We require the statistics of material properties to be statistically homogeneous and mean ergodic. The latter property means that the spatial average equals the ensemble average(2.26)where denotes the ensemble average and is the probability function of the sample *ω*. Note, of course, that is not given by the above equality.

As a particular case, one may assume the statistics of microstructure to be isotropic, which would yield(2.27)i.e. the isotropy of the effective response on the RVE level.

In the following, we partition a square-shaped into 2^{m} subdomains (in two-dimensional, *m*=2; in three-dimensional, *m*=3), *s*=1, 2, …, 2^{m}. Thus, each has the mesoscale . Clearly,(2.28)Henceforth, for simplicity, we work in two dimensions as shown by figure 1*b*; the procedure is totally analogous in three dimensions.

We now define two types of boundary conditions here: namely, the ‘unrestricted’ and ‘restricted’ boundary conditions for and all four 's, respectively.

First, let us consider an ‘unrestricted’ Dirichlet boundary condition: the prescribed pressure on the boundary of :(2.29)In view of the Hill condition, the potential energy can be expressed as(2.30)where is the permeability tensor first appearing in (2.24).

On the other hand, we can also set up a ‘restricted’ Dirichlet boundary condition: the prescribed pressure on the boundary of all four in (2.28):(2.31)which, given the Hill condition, leads to(2.32)Here, is different from of (2.30), whereby ‘r’ denotes the restriction.

We now observe that the solution under the restricted condition (2.31) on all four subdomains is an admissible solution with respect to the unrestricted condition (2.29) on , but not vice versa. Thus, in view of the minimum potential energy principle, we have the inequality: , where are the energies of particular subdomains in (2.28). Henceforth, given (2.28), this inequality relationship becomes(2.33)where pertain to particular subdomains in (2.28). Note that here. Consequently, given the medium's statistical homogeneity and ergodicity, the following relationship is obtained upon taking the ensemble average over the space:(2.34)The inequality relation between tensors in (2.30) means that, when for any , then . By induction, (2.34) leads to a scale-dependent hierarchy (upper bound) of the permeability of porous media flow(2.35)Note that pertains to the smallest scale where Darcy's law still applies; by analogy to elasticity of random media, that tensor may be called a Voigt bound.

Now, let us consider an ‘unrestricted’ Neumann boundary condition: the prescribed velocity on the boundary of ,(2.36)and a ‘restricted’ Neumann boundary condition: the prescribed velocity on the boundary of all four in (2.28),(2.37)Next, consider the energy dissipation rate under (2.36) and (2.37), respectively, via the analogy to complementary energy. Note that the Hill condition is applicable in both cases, and, corresponding to (2.36) and (2.37), the complementary energies are(2.38)Here, is the resistance tensor first appearing in (2.25).

Now, in view of the minimum complementary energy principle (of §2*a*), the solution under the restricted boundary condition (2.37) is admissible with respect to the unrestricted boundary condition (2.36), but not vice versa. Hence,(2.39)where pertain to particular subdomains in (2.28). As before, upon taking the ensemble average over the space of the assumed random porous media, we arrive at the inequality for the resistance tensor on two mesoscales, *δ* and ,(2.40)Clearly, this indicates the upper bound character of resistance tensors with respect to . Since *δ* is arbitrary, from (2.40), we obtain a scale-dependent hierarchy (of lower bounds) of the inverse of the resistance tensor(2.41)Again by analogy to elasticity of random media, tensor pertains to the smallest scale where Darcy's law still applies; and it may be called a Reuss bound.

Noting (2.8), the hierarchies (2.35) and (2.41) can be combined into a single, two-sided hierarchy:(2.42)This relation states that for random heterogeneous porous media, which are statistically homogeneous and ergodic:

the effective permeability tensor is bounded from above by the permeability tensor obtained under the Dirichlet boundary condition, and from below by the permeability tensor obtained under the Neumann boundary condition;

the hierarchy (2.42) shows the scale dependence of apparent permeability tensors under Dirichlet or Neumann boundary condition, which with the increasing mesoscale

*δ*tends to converge towards the effective permeability.

Note that, just like in the case of elastic materials, the hierarchy (2.42) may be extended to arbitrary scales , not only those satisfying (Ostoja-Starzewski 1999).

## 3. The RVE size via computational fluid mechanics

### (a) Direct simulation

The foregoing discussion also indicates a procedure to quantitatively homogenize the heterogeneous medium and determine the RVE size. In the following, we conduct numerical simulations to quantitatively demonstrate the convergent tendency of the hierarchy (2.42). Starting from the energy conservation law in (2.5), we solve Stokes' equation (2.1) numerically in the void of a random model for a given boundary condition, either Dirichlet or Neumann. In order to deal with the geometric complexity of random porous media, we employ a control volume finite-difference method, and to avoid the enhancement of incompressibility by Stokes flow, which causes an unrealistic oscillation of the pressure field in computation, a so-called semi-implicit method for pressure-linked equations revised with collocated equal-order discretization technique is used (Patankar 1980; Chung 2000).

The porous media model investigated in the present research is a two-dimensional flow past a bed of round discs randomly distributed by a planar hard-core Poisson point field (figure 1*a*). The latter condition means that no two discs may touch. In fact, in order to avoid the narrow-neck effects in fluid field, the minimum distance between disc centres is set at 1.1 diameter of the disc. Such effects may be introduced as an extension of our model in future work.

We set up boundary conditions for Darcy's law. First, for the Dirichlet boundary condition with prescribed by (2.20) for a mesoscale specimen at a given *δ*, we must also pose a proper boundary condition for the Stokes flow in the void region. Beside the pressure distribution (2.20), and according to (2.3) and (2.4), we set on the top and bottom of the element, and on the left- and right-hand side of the element for the direct numerical simulation of the void region. This leads to a calculation of the mesoscale permeability tensor .

On the other hand, for the Neumann boundary condition, given in (2.21), also according to (2.3) and (2.4), we set on the top and bottom boundaries of the element. This allows us to calculate the mesoscale resistance tensor .

Figure 2 illustrates the comparison between both kinds of boundary conditions. We can clearly see in figure 2*b* that the pressure distributes nonlinearly on side boundaries in order to maintain the flow within the top and bottom boundary according to the requirement of Neumann boundary condition. Nevertheless, for a sufficiently large mesoscale, the fluctuation of pressure becomes so negligible that the results of the Neumann boundary-value problem begin to coincide with those of the Dirichlet boundary-value problem (except for the obvious non-uniqueness of the Neumann problem). This is the RVE limit we want to assess.

### (b) Numerical results and discussion

To assess the hierarchy (2.42), we proceed with the following steps:

choose a specific mesoscale

*δ*and a nominal porosity;generate a mesoscale realization of the disordered medium, as described in §3

*a*;separately solve the Dirichlet and Neumann boundary-value problems for each , and store the apparent permeability and resistance tensors: and ;

proceed from (i) to (iii) in a Monte Carlo sense over the sample space to generate an ensemble of results at a fixed mesoscale and ;

change the mesoscale and repeat steps (i)–(iv) to find the first moments and ;

change the nominal porosity and repeat steps (i)–(v).

In this procedure, we use *δ*=4, 8, 16, 32, 64 and 80, and nominal porosities 50, 60, 70 and 80%. The corresponding numbers of Monte Carlo realizations are, respectively, 160, 80, 50, 36, 25 and 16.

Figure 3 displays the simulation results of the random medium at 60% porosity in terms of the ensemble average half-traces of the permeability and resistance tensors. The bounding character of these tensors and their convergence to the RVE (effective medium) with the increasing mesoscale are clearly seen. In fact, at *δ*=80, the ensemble averages of both tensors are approximately equal within 0.4%.

Another viewpoint of that convergence to RVE may be offered in terms of the departure of one-half of the ensemble-averaged scalar product of and from unity, for we must have at ,(3.1)The statistical isotropy of the underlying Poisson point field with exclusion dictates the isotropy of not only the effective (RVE level) properties, but also the isotropy of and . Thus, we can work with(3.2)to characterize the above-mentioned departure from the RVE properties; here 1*j* are the indices of particular components of the and tensors. The said departure is shown in figure 4, both for the preceding case of porosity, as well as for the entire range: 50–80%. Clearly, the scaling effect of permeability in porous media depends more strongly on the size of micro-channels than on the disc diameter *d* itself. The rate of convergence to RVE decreases with increasing porosity. For example, at *δ*=32, the RVE size can be considered to be attained at porosity 50%. However, for a higher porosity, say 80%, the RVE size is still not achieved even at *δ*=80.

Of course, the rate of convergence to RVE also depends on the acceptable error. Therefore, in figure 5, we show *δ* with respect to the scale change and porosities. We can see that, for a 10% error, in the medium with 80% porosity, the RVE can be reached at the mesoscale *δ*=56. However, for a 1% error, it cannot be reached even at .

As a final issue, we assess the rate of convergence to the isotropic behaviour. That property is measured using the second invariant of either or , i.e. the radius of the Mohr circle of the second-rank tensor,(3.3)Also, the variance of *r* converges to zero as well. Figure 6 shows the coefficient of variation of *r* (i.e. its standard deviation divided by the mean) in function of the mesoscale. We see that the variance converges to zero faster than the ensemble average in the present study.

## 4. Conclusion

We summarize this study as follows.

We set up the mathematical basis for finding effective permeability and RVE size for porous media with spatially homogeneous and ergodic statistics. By analogy to the existing methodology in thermomechanics of solid random media, we first formulate a Hill–Mandel condition for the Darcy flow velocity and pressure gradient fields. This dictates uniform Neumann and Dirichlet boundary conditions, which, with the help of the minimum potential and complementary energy principles for the Stokes flow, lead to scale-dependent hierarchies on effective (RVE level) permeability.

The scale-dependent hierarchy of mesoscale bounds on the effective (RVE level) permeability tensor is derived: the ensemble average of the permeability tensor is bounded from above by Dirichlet boundary-value problems, and from below by Neumann boundary-value problems. The upper and lower bounds converge towards the RVE response with the window scale increasing.

The hierarchy of bounds is computed for a range of volume fractions of random media made up of planar discs, where the disc centres are generated by a planar hard-core Poisson point field. Overall, it turns out that the higher is the density of random discs—or, equivalently, the narrower are the micro-channels in the system—the smaller is the size of the RVE pertaining to the Darcy law. In particular, the condition in the separation of scales (1.1) is not always necessary, and it is possible that the RVE is attained at .

This paper sets the stage for computation of the RVE size in more general settings: permeability of non-Newtonian fluids; permeability in three-dimensional microstructures; poroelasticity, where the fluid flows in a porous elastic microstructure.

## Acknowledgments

This work has been made possible through the support by the Canada Research Chairs program and the NSERC.

## Footnotes

- Received October 18, 2005.
- Accepted February 27, 2006.

- © 2006 The Royal Society