V.C.1 Permeability and Darcy’s Law
The permeability is the most important physical property
of a porous medium in much the same way as the porosity
is its most important geometrical property.
Some authors define porous media as media with
a nonvanishing permeability [2].
Permeability measures quantitatively the ability
of a porous medium to conduct fluid flow.
The permeability tensor K relates the
macroscopic flow density v¯ to the applied pressure
gradient ∇P or external field F through
where η is the dynamic viscosity of the fluid.
v¯ is the flow rate per unit area of cross section.
Equation (5.54) is known as Darcy’s law.
The permeability has dimensions of an area, and it is
measured in units of Darcy (d).
If the pressure is measured in physical atmospheres
one has 1d=0.9869μm2 while 1d=1.0197μm2
if the pressure is measured in technical atmospheres.
To within practical measuring accuracy one may often
assume 1d=10-12m2.
An important question arising from the fact that K
is dimensionally an area concerns the interpretation of
this area or length scale in terms of the underlying geometry.
This fundamental question has recently found renewed
interest [318, 43, 319, 320, 170, 172, 4].
Unfortunately most answers proposed in these discussions
[319, 318, 320, 4]
give a dynamical rather than geometrical interpretation of
this length scale.
The traditional answer to this basic problem is provided
by hydraulic radius theory [3, 2].
It gives a geometrical interpretation which is based on the
capillary models of section III.B.1, and it will
be discussed in the next section.
The permeability does not appear in the microscopic Stokes
or Navier-Stokes equations.
Darcy’s law and with it the permeability concept can be
derived from microscopic Stokes flow equations using
homogenization techniques
[268, 269, 270, 38, 321, 271]
which are asymptotic expansions in the ratio of
microscopic to macroscopic length scales.
The derivation will be given in section V.C.3
below.
The linear Darcy law holds for flows at low Reynolds numbers
in which the driving forces are small and balanced only by the
viscous forces.
Various nonlinear generalizations of Darcy’s law have also
been derived using homogenization or volume averaging methods
[268, 1, 269, 322, 321, 38, 271, 323, 324, 325].
If a nonlinear Darcy law governs the flow in a given experiment
this would appear in the measurement as if the permeability
becomes velocity dependent.
The linear Darcy law breaks down also if the flow becomes too slow.
In this case interactions between the fluid and the pore walls
become important.
Examples occur during the slow movement of polar liquids or electrolytes
in finely porous materials with high specific internal surface.
V.C.2 Hydraulic Radius Theory
The hydraulic radius theory or Carman-Kozeny model
is based on the geometrical models of capillary
tubes discussed above in section III.B.1.
In such capillary models the permeability can be
obtained exactly from the solution of
the Navier-Stokes equation (4.9)
in the capillary.
Consider a cylindrical capillary tube of length L
and radius a directed along the x-direction.
The velocity field vr for creeping laminar
flow is of the form vr=vrex
where ex denotes a unit vector along the pipe, and
r measures the distance from the center of the pipe.
The pressure has the form Pr=Pxex.
Assuming “no slip” boundary conditions, va=0, at the tube
walls one obtains for vr the familiar Hagen-Poiseuille
result [326]
Px | = | P0-P0-PLxL |
| (5.55) |
vr | = | P0-PL4ηLa2-r2 |
| (5.56) |
with a parabolic velocity and linear pressure profile.
The volume flow rate Q is obtained through integration as
Q=∫0avr2πrdr=πa48ηP0-PLL. |
| (5.57) |
Consider now the capillary tube model of section III.B.1
with a cubic sample space 𝕊 of sidelength L.
The pore space ℙ consists of N nonintersecting capillary
tubes of radii ai and lengths Li distributed according
to a joint probability density Πa,L.
The pressure drop must then be calculated over the length Li
and thus the right hand side of (5.57) is
multiplied by a factor L/Li.
Because the tubes are nonintersecting the volume flow Qi
through each of the tubes can be added to give the macroscopic
volume flow rate per unit area v¯=1/L2∑i=1NQi.
Thus the permeability of the capillary tube model
is simply additive, and it reads
k=π8L∑i=1Nai4Li. |
| (5.58) |
Dimensional analysis of (5.58), (3.58)
and (3.59) shows that kS2/ϕ3 is
dimensionless.
Averaging (5.58) as well as (3.58)
and (3.59) for the porosity and specific internal
surface of the capillary tube model yields the relation
where the mixed moment ratio
C=L2a4LaL2a2L3 |
| (5.60) |
is a dimensionless number, and the angular brackets denote
as usual the average with respect to Πa,L.
The hydraulic radius theory or
Carman-Kozeny model is obtained from a mean field
approximation which assumes fx≈fx .
The approximation becomes exact if the distribution is sharply
peaked or if Li=L and ai=a for all N.
With this approximation
the average permeability k may be rewritten
in terms of the average hydraulic radius RH defined in
(3.66) as
k≈ϕ2𝒯2ϕ2S2≈ϕ2𝒯2ϕ2S2≈ϕRH22𝒯2≈ϕa28𝒯2 |
| (5.61) |
where 𝒯=L/L is the average of the tortuosity
defined above in (3.62).
Equation (5.61) is one of the main results of
hydraulic radius theory.
The permeability is expressed as the square of an
average hydraulic radius RH, which is related to the
average “pore width” as RH=a/2.
It must be stressed that hydraulic radius theory is not exact
even for the simple capillary tube model because in general
RH≠ϕ/S and C≠1/𝒯2.
However, interesting exact relations for the average permeability
can be obtained from (5.59) and (5.60)
in various special cases without employing the
mean field approximation of hydraulic radius theory.
If the tube radii and lengths are independent then the distribution
factorizes as Πa,L=ΠaaΠLL.
In this case the permeability may be written as
k=121/𝒯𝒯a4a2a23ϕ3S2=ϕ81/𝒯𝒯a4a2 |
| (5.62) |
where 𝒯 is the average of the tortuosity factor
defined in (3.62).
The last equality interprets k in terms of the
microscopic effective cross section a4/a2
determined by the variance and curtosis of the distribution
of tube radii.
Further specialization to the cases Li=l or ai=a is
readily carried out from these results.
Finally it is of interest to consider also the capillary slit
model of section III.B.1.
The model assumes again a cubic sample of side length L
containing a pore space consisting of parallel slits with
random widths governed by a probability density Πb.
For flat planes without undulations the analogue of tortuosity
is absent.
The average permeability is obtained in this case as
which has the same form as (5.59) with a constant
C=b3/b3.
The prefactor 1/3 is due to the different shape of the
capillaries, which are planes rather than tubes.
V.C.3 Derivation of Darcy’s Law from Stokes Equation
The previous section has shown that Darcy’s law arises in the
capillary models.
This raises the question whether it can be derived more generally.
The present section shows that Darcy’s law can be obtained from
Stokes equation for a slow flow.
It arises to lowest order in an asymptotic expansion
whose small parameter is the ratio of microscopic to macroscopic
length scales.
Consider the stationary and creeping (low Reynolds number)
flow of a Newtonian incompressible fluid through a porous
medium whose matrix is assumed to be rigid.
The microscopic flow through the pore space ℙ is governed
by the stationary Stokes equations for the velocity vr
and pressure Pr
ηΔvr+F-∇Pr | = | 0 |
| (5.64) |
∇T⋅vr | = | 0 |
| (5.65) |
inside the pore space, ℙ∋r, with no slip boundary condition
for r∈∂ℙ.
The body force F and the dynamic viscosity η are
assumed to be constant.
The derivation of Darcy’s law assumes that the pore space ℙ
has a characteristic length scale l which is small compared to
some macroscopic scale L.
The microscopic scale l could be the diameter of grains,
the macroscale L could be the diameter of the sample 𝕊
or some other macroscopic length such as the diameter of
a measurement cell or the wavelength of a seismic wave.
The small ratio ε=l/L provides a small parameter for
an asymptotic expansion.
The expansion is constructed by assuming that all properties
and fields can be written as functions of two new space variables
x,y which are related to the original space variable
r as x=r and y=r/ε.
All functions fr are now replaced with functions fx,y
and the slowly varying variable x is allowed to vary
independently of the rapidly varying variable y.
This requires to replace the gradient according to
∇fr=∇fr,r/ε=∇xfx,y+1ε∇yfx,y |
| (5.67) |
and the Laplacian is replaced similarly.
The velocity and pressure are now expanded in ε where
the leading orders are chosen such that the solution is not
reduced to the trivial zero solution and the problem remains
physically meaningful.
In the present case this leads to the expansions [268, 280, 271]
vr | = | ε2v0x,y+ε3v1x,y+… |
| (5.68) |
Pr | = | P0x,y+εP1x,y+… |
| (5.69) |
where x=r and y=r/ε.
Inserting into (5.64), (5.65) and (5.66)
yields to lowest order in ε the system of equations
∇yP0x,y | = | 0in ℙ |
| (5.70) |
∇yT⋅v0 | = | 0in ℙ |
| (5.71) |
ηΔyv0-∇yP1-∇xP0+F | = | 0in ℙ |
| (5.72) |
∇xT⋅v0+∇yT⋅v1 | = | 0in ℙ |
| (5.73) |
v0 | = | 0on ∂ℙ |
| (5.74) |
in the fast variable y.
It follows from the first equation that P0x,y depends
only on the slow variable x, and thus it appears as an additional
external force for the determination of the dependence of
v0x,y on y from the remaining equations.
Because the equations are linear the solution v0x,y
has the form
v0x,y=∑i=13Fi-∂P0∂xiuix,y |
| (5.75) |
where the three vectors uix,y (and the scalars Qix,y)
are the solutions of the three systems (i=1,2,3)
∇yT⋅ui | = | 0in ℙ |
| (5.76) |
ηΔyui-∇yQi-eyi | = | 0in ℙ |
| (5.77) |
ui | = | 0on ∂ℙ |
| (5.78) |
and eyi is a unit vector in the direction of the
yi-axis.
It is now possible to average v0 over the fast variable y.
The spatial average over a convex set 𝕂 is defined as
v¯0x;𝕂=1V𝕂∫v0x,yχ𝕂x,yd3y |
| (5.79) |
where 𝕂 is centered at x and
χ𝕂x,y=χ𝕂r,r/ε equals 1 or 0 depending
upon whether r∈𝕂 or not.
The dependence on the averaging region 𝕂 has been indicated
explicitly.
Using the notation of (2.20) the average over all
space is obtained as the limit lims→∞v¯0x;s𝕂=v¯0x.
The function P0 need not to be averaged as it depends only on the slow
variable x.
If v0 is constant then v0¯x=v0ϕ¯x which
is known as the law of Dupuit-Forchheimer [1].
Averaging (5.75) gives Darcy’s law (5.54)
in the form
v¯0x;𝕂=Kx;𝕂ηF-∇xP0x |
| (5.80) |
where the components kijx;𝕂=Kx;𝕂ij of
the permeability tensor K are expressed in terms of
the solutions ujx;𝕂 to (5.76)–(5.78)
within the region 𝕂 as
Kx;𝕂ij=uj¯x;𝕂i. |
| (5.81) |
The permeability tensor is symmetric and positive definite
[268].
Its dependence on the configuration of the pore space ℙ
and the averaging region 𝕂 have been made explicit because
they will play an important role below.
For isotropic and strictly periodic or stationary media the
permeability tensor reduces to a constant independent of x.
For (quasi-)periodic microgeometries or (quasi-)stationary random
media averaging eq. (5.73) leads to the additional
macroscopic relation
Equations (5.80) and (5.82)
are the macroscopic laws governing the microscopic
Stokes flow obeying (5.64)–(5.66)
to leading order in ε=l/L.
The importance of the homogenization technique illustrated
here in a simple example lies in the fact that it provides
a systematic method to obtain the reference problem for an
effective medium treatment.
Many of the examples for transport and relaxation in
porous media listed in chapter IV
can be homogenized using a similar technique [268].
The heterogeneous elliptic equation (4.2)
is of particular interest.
The linear Darcy flow derived in this section can be cast into
the form of (4.2) for the pressure field.
The permeability tensor may still depend on the slow variable x,
and it is therefore of interest to iterate the homogenization
procedure in order to see whether Darcy’s law becomes again
modified on larger scales.
This question is discussed next.
V.C.4 Iterated Homogenization
The permeability Kx for the macroscopic Darcy flow
was obtained from homogenizing the Stokes equation by
averaging the fast variable y over a region 𝕂.
The dependence on the slow variable x allows for
macroscopic inhomogeneities of the permeability.
This raises the question whether the homogenization
may be repeated to arrive at an averaged description
for a much larger megascopic scale.
If (5.80) is inserted into (5.82)
and F=0 is assumed the equation for the macroscopic pressure
field becomes
which is identical with (4.2).
The equation must be supplemented with boundary
conditions which can be obtained from the requirements
of mass and momentum conservation at the boundary of the
region for which (5.83) was derived.
If the boundary marks a transition to a region with
different permeability the boundary conditions require
continuity of pressure and normal component of the
velocity.
Equation (5.83) holds at length scales L
much larger than the pore scale l, and much larger than
diameter of the averaging region 𝕂.
To homogenize it one must therefore consider length scales
ℒ much larger than l such that
is fulfilled.
The ratio δ=L/ℒ is then a small parameter in terms of
which the homogenization procedure of the previous section can
be iterated.
The pressure is expanded in terms of δ as
Px=P0s,z+δP1s,z+… |
| (5.85) |
where now s=x is the slow variable, and z=s/δ
is the rapidly varying variable.
Assuming that the medium is stationary, i.e. that Kz
does not depend on the slow variable s,
the result becomes [268, 280, 271]
where P0s is the first term in the expansion
of the pressure which is independent of z, and the
tensor K¯ has components
K¯ij=kijz+∑l=13kilz∂Qjz∂zl¯ |
| (5.87) |
given in terms of three scalar fields Qj(j=1,2,3) which
are obtained from solving an equation of the form
-∑i,j∂∂zikijz∂Qkz∂zj=∑i∂kikz∂zi |
| (5.88) |
analogous to (5.76)–(5.78) in the homogenization
of Stokes equation.
If the assumption of strict stationarity is relaxed the
averaged permeability depends in general on the slow variable,
and the homogenized equation (5.86) has then the
same form as the original equation (5.83).
This shows that the form of the macroscopic equation does not
change under further averaging.
This highlights the importance of the averaged permeability
as a key element of every macroscopically homogeneous description.
Note however that the averaged tensor K¯
may have a different symmetry than the original permeability.
If Kx=kx1 is isotropic (1 denotes
the unit matrix) then K¯ may become anisotropic
because of the second term appearing in (5.87).
V.C.5 Network Model
Consider a porous medium described by equation
(5.83) for Darcy flow with a stationary
and isotropic local permeability function
Kx=kx1.
The expressions (5.87) and (5.88) for the
the effective permeability tensor K¯ are difficult to
use for general random microstructures.
Therefore it remains necessary
to follow the strategy outlined in section V.A.2
and to discretize (5.83) using a finite difference
scheme with lattice constant L.
As before it is assumed that l≪L≪ℒ where l is
the pore scale and ℒ is the system size.
The discretization results in the linear network equations (5.5)
for a regular lattice with lattice constant L.
To make further progress it is necessary to specify
the local permeabilities.
A microscopic network model of tubes results from choosing
the expression
for a cylindrical capillary tube of radius a
and length ℓ in a region of size L.
The parameters a and ℓ must obey the geometrical
conditions a≤L/2 and ℓ≥L.
In the resulting network model each bond
represents a winding tube with circular cross
section whose diameter and length fluctuate
from bond to bond.
The network model is completely specified by assuming
that the local geometries specified by a and ℓ
are independent and identically distributed random
variables with joint probability density Πa,ℓ.
Note that the probability density Πa,ℓ depends
also on the discretization length through the constraints
a≤L/2 and ℓ≥L.
Using the effective medium approximation to the network
equations the effective permeability k¯ for this
network model is the solution of the selfconsistency equation
∫L∞∫0L/2πa4-8Lℓk¯πa4+16Lℓk¯Πa,ℓdadℓ=0 |
| (5.90) |
where the restrictions on a and ℓ are reflected
in the limits of integration.
In simple cases, as for binary or uniform distributions,
this equation can be solved analytically, in other
cases it is solved numerically.
The effective medium prediction agrees well with
an exact solution of the network equations
[231].
The behaviour of the effective permeability depends
qualitatively on the fraction p of conducting tubes
defined as
p=1-limε→0∫0εΠada |
| (5.91) |
where Πa=∫L∞Πa,ℓdℓ.
For p>1/3 the permeability is positive while
for p<1/3 it vanishes.
At p=pc=1/3 the network has a percolation transition.
Note that p≠ϕ¯ is not related to the
average porosity.
V.C.6 Local Porosity Theory
Consider, as in the previous section, a porous medium described
by equation (5.83) for Darcy flow with a stationary
and isotropic local permeability function Kx=kx1.
A glance at section III shows that the one cell
local geometry distribution defined in (3.45) are particularly
well adapted to the discretization of (5.83).
As before the discretization employs a cubic lattice with lattice
constant L and cubic measurement cells 𝕂 and yields a
local geometry distribution μϕ,S;𝕂.
It is then natural to use the Carman equation (5.59)
locally because it is often an accurate description
as illustrated in Figure 23.
The straight line in Figure 23 corresponds
to equation (5.59).
The local percolation probabilities defined in section
III.A.5.d complete the description.
Each local geometry is characterized by its local
porosity, specific internal surface and a binary random
variable indicating whether the geometry is percolating
or not.
The selfconsistent effective medium equation now reads
∫0∞∫013Cϕ3λϕ,S;𝕂μϕ,S;𝕂Cϕ3+4S2k¯dϕdS=1 |
| (5.92) |
for the effective permeability k¯.
The control parameter for the underlying percolation transition
was given in (3.47) as
pL=∫0∞∫01λϕ,S;𝕂μϕ,S;𝕂dϕdS |
| (5.93) |
and it gives the total fraction of percolating local geometries.
If the quantity
k0=∫0∞∫012S2Cϕ3λϕ,S;𝕂μϕ,S;𝕂dϕdS-1 |
| (5.94) |
is finite then the solution to (5.92) is given approximately as
for p>pc=1/3 and as k¯=0 for p<pc.
This result is analogous to (5.48) for the
electrical conductivity.
Note that the control parameter for the underlying percolation
transition differs from the bulk porosity p≠ϕ¯.
To study the implications of (5.92) it is necessary
to supply explicit expressions for the local geometry distribution
μϕ,S;L.
Such an expression is provided by the local porosity reduction
model reviewed in section III.B.6.
Writing the effective medium approximation for the number n¯
defined in (3.87) and using equations (3.86)
and (3.88) it has been shown that the effective
permeability may be written approximately as [170]
where the exponent β depends on the porosity reduction
factor r and the type of consolidation model characterized by
(3.88) as
If all local geometries are percolating, i.e. if λ=1,
then the effective permeability depends algebraically on the
bulk porosity ϕ¯ with a strongly nonuniversal
exponent β.
This dependence will be modified if the local percolation
probabilityλϕ¯ is not constant.
The large variability is consistent with experience from
measuring permeabilities in experiment.
Figure 24 demonstrates the large data scatter
seen in experimental results.
While in general small permeabilities correlate with small
porosities the correlation is not very pronounced.