PHYSICAL REVIEW E VOLUME 54, NUMBER 1 JULY 1996

Tortuous flow in porous media

A. Koponen, M. Kataja, and J. Timonen
Department of Physics, University of Jyväskylä, PL 35, SF-40351 Jyväskylä, Finland
~Received 2 January 1996!
The concept of tortuosity of fluid flow in porous media is discussed. A lattice-gas cellular automaton method
is applied to solve the flow of a Newtonian uncompressible fluid in a two-dimensional porous substance
constructed by randomly placed rectangles of equal size and with unrestricted overlap. A clear correlation
between the average tortuosity of the flow paths and the porosity of the substance has been found. @S1063-
651X~96!05006-4#

PACS number~s!: 47.15.2x, 47.55.2t

I. INTRODUCTION cy’s law which states that the average flux q of fluid is pro-
portional to the gradient of the phase averaged fluid pressure
A common characteristic of any material transport in po- p,
rous media, such as fluid flow or electric current, is that the
actual path followed by the transported material is micro- k
scopically very complicated, or ‘‘tortuous’’ @1–4# ~‘‘micro- q52 “p. ~1!
m
scopic’’ here means the size scale of the average pore size of
the substance!. The concept of tortuosity is often introduced
in the context of solving the closure problem for transport in Here, m is the dynamic viscosity of the fluid and k is the flow
porous media, i.e., in deriving the macroscopic transport permeability. In this paper, we shall first discuss the possible
equations in terms of averaged quantities alone. A usual definitions of tortuosity which appear in deriving Darcy’s
method of deducing, e.g., the appropriate form of the drag law in the framework of capillary models. The definition is
force between fluid and the solid matrix, is to use some sim- then generalized to random porous media. We then use a
plified model of the porous material, such as the capillary lattice-gas cellular automaton method @6–8# to find a numeri-
model, and to generalize the results for more realistic mate- cal correlation between the tortuosity and the porosity for a
rials. This generalization may be attempted by introducing an two-dimensional porous material which consists of randomly
additional parameter which is supposed to take care of the positioned freely overlapping rectangles. The advantage of
more complicated transport paths neglected in the model. In using the lattice-gas methods for this purpose is their geo-
fact tortuosity is an example of such a parameter. As a physi- metric versatility, which makes them very useful in simulat-
cal quantity, it can be defined in various ways. Perhaps the ing flows in irregular geometries @9–11#. The results ob-
most intuitive and straightforward definition is that of the tained can be applied, e.g., in inferring correlations between
ratio of the average length of true flow paths to the length of the permeability coefficient and the relevant macroscopic
the system in the direction of the macroscopic flux. Notice quantities that characterize the porous substance @12#. Find-
that by this definition, tortuosity depends not only on the ing such correlations is especially important in the case of
microscopic geometry of the pores, but also on the transport flow through soft porous materials when flow can induce a
mechanism under consideration. strain to the solid matrix, and thereby locally affect the value
Tortuosity could also be defined without reference to a of the permeability coefficient k @13#.
specific transport mechanism. This could be done, for ex-
ample, by considering the shortest continuous paths between II. TORTUOSITY OF FLOW IN POROUS MEDIA
any two points within the pore space @5#. The advantage of
this definition is that the tortuosity parameter thus defined Darcy’s law Eq. ~1! can easily be derived within the
will exclusively characterize the porous substance itself. simple capillary theory by Kozeny, in which the porous me-
When considering tortuosity in the context of transport phe- dium is envisaged as a layer of solid material with straight
nomena, it seems quite more natural, however, to utilize the parallel tubes of a fixed cross-sectional shape intersecting the
flux associated with the actual transport mechanism in the sample. Within this model, the permeability is explicitly
definition of tortuosity. given as k5 f 3 /cs 2 , where f is the porosity, s is the specific
Moreover, it is possible to define tortuosity even without a surface area, i.e., the pore surface area per unit volume of
direct reference to the lengths of the transport paths by con- porous material, and c is a structural parameter that depends
sidering the local deviations in the direction of the micro- on the cross section of the capillaries, for cylindrical capil-
scopic flux from the direction of the mean flux. ~This ap- laries c52. The simplest way to introduce tortuosity in the
proach will be discussed in some detail below.! capillary model is to allow the tubes to be inclined in such a
In what follows, we shall concentrate on the concept of way that the axes of the capillaries form a fixed angle u with
tortuosity associated with the flow of a Newtonian fluid the normal of the surface of the material ~while the azimuthal
through a random porous medium at a low Reynolds num- angle of the tubes is randomly distributed!. In this case per-
ber. In macroscopic terms, such a flow is governed by Dar- meability becomes

1063-651X/96/54~1!/406~5!/$10.00 54 406 © 1996 The American Physical Society

54 TORTUOUS FLOW IN POROUS MEDIA 407

f3 Notice first that, within the capillary model, the product of

k5 , ~2! A i , the area of the intersection of the ith tube with plane A,
c t 2s 2
and the average flow velocity in the ith capillary v i is inde-
where the tortuosity factor t 51/cosu can be given in terms pendent of i ( v i A i 5 p R 4 u “p u /8m ). Using this result we find
of the tube length L e and the thickness of the medium L as that
N
t 5L e /L. ~3!
1
N ( t̃ i v i A i * A t̃ v dA
„Some authors prefer to define tortuosity as t 5(L e /L) or as
( t̃ i 5
2 i51
5 . ~6!
the inverses of these two definitions @3,4#. In this paper, we N N
* A v dA
( v iA i
i51
shall use definitions analogous to Eq. ~3!. Thus for the tor- i51
tuosity defined here t>1.…
For flow in random porous media, one can replace the Here, v 5uv~r!u is the tangential velocity of the fluid at point
‘‘tube length’’ L e by the average length of the flow paths of r, and t̃ 5 t̃ ~r! is the ratio of the length of the flow line
a fluid particle through the sample. At least two possible passing through the point r to the thickness of the sample
alternatives for taking this average can be considered @3#. @t̃ ~r! and v are defined to be zero inside the solid phase#. A
One may average over the actual lengths of the flow lines similar result is valid for the sum containing 1/t̃ i . This sug-
themselves, disregarding thereby the fact that fluid particles gests the following definitions,
move along these flow lines at different velocities. Another
way of averaging is over the lengths of the flow lines of all * A t̃ v dA * V t̃ v dV
t 15 5 , ~7!
fluid particles passing through a given cross section during a * A v dA * V v dV
given period of time. This leads to flux weighted averaging.
The first alternative is suitable at least for pistonlike flows, 1 1
such as molecular diffusion and electric current @3#. The lat- *A v dA * V v dV
t̃ t̃
ter alternative appears more natural when fluid flow in po- 1/t 21 5 5 , ~8!
rous media is considered. * A v dA * V v dV
In order to gain more insight into the definition ~or defi- where V is the volume of the porous sample. The latter forms
nitions! of tortuosity of flow in porous media, we shall also of Eqs. ~7! and ~8! follow from the fact that surface integrals
consider a solid material of thickness L, intersected by N do not depend on the position ~x coordinate! of surface A
cylindrical capillaries per unit transverse area. We assume inside the sample. The tortuosity factor as determined from
the capillaries are straight and of equal radius R, but allow Eq. ~7! can be interpreted as the average of the relative
for a randomly varying angle between them and the x axis, lengths of the flow lines of all fluid elements ~with a fixed
which is perpendicular to the surfaces of the material. For the volume! passing through a given cross section during a given
ith capillary of length L i we define, in accordance with Eq. period of time. The latter definition Eq. ~8! corresponds to
~3!, t̃ i 5L i /L. Next, flow through the capillaries is induced the average of inverse lengths of the same flow lines. The
by applying a pressure difference Dp across the sample. tortuosity factor ~5! appearing in Darcy’s law Eq. ~4! can
Solving the Navier-Stokes equation for each capillary sepa- now be expressed in a generalized form which is applicable
rately, we can determine the average flux through the in random porous materials,
sample, with the result
N t 2 5 t 1 t 21 . ~9!
1
f “p N i51
3
1/t̃ i ( Equations ~7! through ~9! do not, however, provide the only
q52 2 , ~4! way of generalizing the results of the capillary model to the
2s m 1 N
t̃
N i51 i ( random media. For example, in the case of capillary systems,
the tortuosity factor t21 of Eq. ~8! is in fact equal to the ratio

where “ p5(Dp/L)êx is the phase averaged pressure gradi- ^ u vu &

t v[ , ~10!
ent. ^ v x&
Comparison of Eq. ~4! with Eqs. ~1! and ~2! suggests a
definition for the tortuosity within this capillary model in the where uvu is the absolute value of the local flow velocity, v x
form is the x component of that velocity, and ^ & denotes the spa-
tial average over the pore space. Notice that Eq. ~10! is remi-
N
1 niscent of the hypothesis made by Carman in Ref. @1# that
N ( t̃ i
i51
L e /L5V̄/ū e , where V̄ is the average tangential velocity in a
t 5
2
N . ~5! tortuous capillary, L e is the length of that capillary, ū e is the
1
( 1/t̃ i
mean value of the projection of flow velocity on the straight
N i51 line connecting the two ends of the capillary, and L is the
length of that line. According to Eq. ~10!, t v is solely deter-
This definition can be expressed in a form that is more suit- mined by fluctuations of the local flow field around the di-
able for generalization for random porous media by convert- rection of the average flux, and has no direct connection with
ing the sums into integrals over an arbitrary plane A which is the length of the actual flow paths. In deriving the above
perpendicular to the x axis ~direction of the average flux!. results, we have assumed that the radius of the capil-
408 A. KOPONEN, M. KATAJA, AND J. TIMONEN 54

laries is fixed while their lengths may vary. The results are,
however, valid also in the case of varying capillary radii,
provided that these and the lengths of the capillaries are un-
correlated.
Yet another possibility, which may be encountered in
other kinds of models, is to define the tortuosity as an aver-
age of the lengths of flow lines squared @3,4#. Analogously to
Eqs. ~7! and ~8!, we may then define

* V v t̃ 2 dV
~ t 2 !25 , ~11!
* V v dV

1
* Vv dV
1 t̃ 2
5 . ~12!
~ t 22 ! 2 * V v dV

At this point, we shall not try to select a preferred defini-

tion of tortuosity among the ones discussed above. Instead,
FIG. 1. Flow lines through a two-dimensional random porous
we will use the lattice-gas word-cellular automaton method
medium.
to find numerical correlations between these tortuosity fac-
tors and the porosity of a two-dimensional random porous N
medium. In a forthcoming publication, these results will be
used to find correlations between the flow permeability and ( t̃ ~ ri !v~ ri !
i51
the macroscopic characteristics of the medium @12#. t 1' N . ~13!
( v~ ri !
i51
III. LATTICE-GAS SIMULATIONS

We solved numerically the two-dimensional flow in a ran- Here, N51000 is the number of flow lines found for each
dom porous medium using the FHP-III lattice-gas model @6# configuration, t̃ ~ri )5L~ri )/L x is the tortuosity of the ith flow
in a discrete triangular mesh of 1003100 lattice sites. The line, L x is the length of the lattice in the x direction, and
two-dimensional porous medium was constructed by random v ~ri ! is the averaged tangential velocity of the fluid at the
positioning of rectangles of an equal size of 10310 lattice starting point. Tortuosities t21, t2 and t22 were determined
sites and with unrestricted overlap. The porosity f of the using expressions similar to Eq. ~13!, while t v was deter-
medium was defined as the ratio of the number of unoccu- mined directly from the averaged velocity field according to
pied sites to the number of all lattice sites. The number of Eq. ~10!.
rectangles K in the lattice varied between 10 and 68, which In Fig. 1 we show an example of the flow line fields for a
corresponds to porosities ranging from 0.9 to 0.5. ~It is configuration of 30 rectangles corresponding to a porosity
straightforward to show that, with the numbers given above, f50.74. The vertical distances between the flow lines of this
the expectation value of the porosity for a given K is figure are determined such that the flux between neighboring
^ f & 50.99K .! The number of fluid particles per lattice site flow lines is constant. The tortuosity t, as determined from
was 3.5 which provides the best approximation for the solu- Eq. ~9!, is t51.2 for this particular configuration. The values
tion of the linearized Navier-Stokes equation ~creeping flow! of tortuosities t1 and t21 @see Eqs. ~7! and ~8!# were found to
within the present method @6#. The fluid was forced to move be very close to each other ~see Fig. 3!. In Fig. 2 we show
in the positive x direction by applying an external force on the calculated tortuosity t 5 At 1 t 21 as a function of porosity
the particles @7#. Periodic boundary conditions were imposed f. The small dots give the values of f and t for the 1080
on the lattice in both directions. individual configurations. The large dots with error bars
Simulations were carried out for about 35 configurations show the mean value and the standard deviation of t and f at
for each K, and the total number of different configurations each value of K. Within the porosity range covered by these
was 1080. A single configuration used about 2.2 hours CPU simulations, the dependence on porosity of tortuosity t is
time on a Dec 3000 workstation. For each configuration the approximately linear. The solid line shown in Fig. 2 is a fit
system was first allowed to saturate for 40 000 time steps by
which was found to yield a stationary flow pattern. The local
velocities of particles were then averaged over 400 000 time t 50.8~ 12 f ! 11. ~14!
steps in order to ensure an undisturbed and smooth flow
velocity field. Flow lines starting from randomly chosen In Fig. 3 we compare the simulated values of the tortuosi-
places ri in the pore space were found by interpolating the ties defined by Eqs. ~7!–~12!. The curves shown in this fig-
time-averaged flow velocity field. The length L~ri ! of each ure are fits to the determined points which are not shown.
flow line ~within one length period in the direction x! was For t v the fit is parabolic while for the other tortuosities it is
then computed and the tortuosity t1 , e.g., was calculated linear. It is evident that the different definitions give approxi-
using the following approximation for the volume integral of mately the same qualitative dependence on porosity, and that
Eq. ~7!, the numerical values do not dramatically differ from each
54 TORTUOUS FLOW IN POROUS MEDIA 409

FIG. 2. The calculated tortuosity t @see Eq. ~9!# as a function of FIG. 3. Comparison of tortuosity factors defined by Eqs. ~7!–
porosity f. Small dots are the results of numerical solutions for ~12!.
1080 individual configurations equivalent to that shown in Fig. 1.
numerical uncertainties was estimated to be below 5%, even
Large symbols with error bars indicate the mean values and stan-
at the lowest porosities where failing of flow lines was most
dard deviations of porosity and tortuosity for each fixed number of
pronounced.
solid blocks.
It is evident that, as a physical quantity, tortuosity is not
uniquely defined. The preferred definition depends on the
other in this porosity range. Notice finally that the presented context and on the model being used. Our simulations sug-
results are in full agreement with an approximate upper limit gest, however, that the model dependence is quite small, at
of tortuosity 1.6 which follows from a model of randomly least for a two-dimensional flow at relatively high porosities.
oriented connected straight tubes @3,4# in two dimensions. Usually the purpose of introducing tortuosity, as a parameter
in macroscopic theories dealing with transport in porous me-
IV. DISCUSSION dia, is to add an additional degree of freedom to account for
the rather complex structure of real porous materials. As
We have used the lattice-gas simulation method for solv- such, tortuosity can hardly be expected to provide more than
ing a low Reynolds number flow in a two-dimensional ma- a qualitative description of the true transport dynamics in
trix formed by randomly placed fully overlapping rectangles. these complex structures. The smallness of the differences
Numerical uncertainties were found to be reasonably small between the numerical values of this quantity, arising from
provided that long enough simulation and averaging times its various plausible definitions, seems to indicate that tortu-
were used to ensure stationary states and smooth velocity osity indeed is a useful concept.
profiles. For a given obstacle configuration the tortuosities The determined interrelation Eq. ~14! of porosity and tor-
calculated with different lattice resolutions were always tuosity can be applied, e.g., in inferring relations between
found to be close to each other, and no systematic resolution permeability and porosity @12#. The basic limitation in doing
effects were seen. this is that the present simulations are two dimensional. For a
In some cases the procedure for finding a flow line pass- three-dimensional flow around nonelongated particles, the
ing a given starting point failed since the flow line hit a solid relation between tortuosity and porosity may not be of the
wall. The contribution from such flow lines was neglected. same form as for two-dimensional objects considered here.
~See, e.g., the second flow line from the top in Fig. 1!. Such Also, two dimensionality restricts the useful configurations
cases were most frequently found in blocked areas where the to those with quite a high porosity. This is due to the perco-
residual fluctuating component of the velocity was relatively lation threshold, which is approximately at f50.33 for ran-
large as compared to the averaged flow velocity. The total domly placed and freely overlapping obstacles ~whose length
flux associated with the failed flow lines and thus their total to width ratio is approximately 1! in two dimensions @14#.
weight in the tortuosity equations was however small. Mak- Close to this porosity, simulations with the present method
ing a conservative assumption that the true lengths of the fail. We therefore expect that the results shown here will be
failed flow lines would differ by at most 30% from those of most directly applicable to flow in fibrous porous media with
the successful flow lines, the error caused by this and other a high porosity.

@1# P. C. Carman, Trans. Inst. Chem. Eng. London 15, 150 ~1937!. @5# P. M. Adler, Porous Media. Geometry and Transports
@2# A. E. Scheidegger, The Physics of Flow Through Porous Me- ~Butterworth-Heinemann, Boston, 1992!.
dia ~Macmillan, New York, 1957!. @6# U. Frisch, D. d’Humiéres, B. Hasslacher, P. Lallemand, Y.
@3# J. Bear, Dynamics of Fluids in Porous Media ~Dover, New Pomeau, and J.-P. Rivet, Complex Syst. 1, 649 ~1987!.
York, 1972!. @7# L. P. Kadanoff, G. R. McNamara, and G. Zanetti, Phys. Rev.
@4# F. A. L. Dullien, Porous Media. Fluid Transport and Pore A 40, 4527 ~1989!.
Structure ~Academic, San Diego, 1979!. @8# D. H. Rothman and S. Zaleski, Rev. Mod. Phys. 66, 1417
410 A. KOPONEN, M. KATAJA, AND J. TIMONEN 54

~1994!. man, and B. J. Travis, Physica D 47, 72 ~1991!.

@9# D. H. Rothman, Geophys. 53, 509 ~1988!. @12# A. Koponen, M. Kataja, and J. Timonen ~unpublished!.
@10# A. Cancelliere, C. Chang, E. Foti, D. H. Rothman, and S. @13# M. Kataja, A. Rybin, and J. Timonen, J. Appl. Phys. 72, 1271
Succi, Phys. Fluids A 2, 2085 ~1990!. ~1992!.
@11# S. Chen, K. Diemer, G. D. Doolen, K. Eggert, C. Fu, S. Gut- @14# G. E. Pike and C. H. Seager, Phys. Rev. B 10, 1421 ~1974!.