Table Of ContentTortuosity–porosity relation in the porous media flow
Maciej Matyka,1 Arzhang Khalili,2,3 and Zbigniew Koza1
1Institute of Theoretical Physics, University of Wrocl aw, pl. M. Borna 9, 50-204 Wrocl aw, Poland
2Max Planck Institute for Marine Microbiology, Celsiusstrasse 1, D-28359 Bremen, Germany
3Jacobs University Bremen, Campus Ring 1, D-28759 Bremen, Germany
(Dated: February 2, 2008)
Westudynumericallythetortuosity–porosityrelationinamicroscopicmodelofaporousmedium
arranged as a collectin of freely overlapping squares. It is demonstrated that the finite-size effects
and the discretization errors, which were ignored in previous studies, may cause significant under-
8 estimation of tortuosity. The simple tortuosity calculation method proposed here eliminates the
0 need for using complicated, weighted averages. The numerical results presented here are in good
0 agreementwithanempiricalrelationbetweentortuosity(T)andporosity(φ)givenbyT−1∝lnφ,
2
that was found by others experimentally in granule packings and sediments. This relation can be
n also written as T −1 ∝ RS/φ with R and S denoting the hydraulic radius of granules and the
a specific surface area, respectively.
J
2 PACSnumbers: 47.56.+r,47.15.G-,91.60.Np
2
] I. INTRODUCTION This effect can be described by a dimensionless parame-
n ter T called hydraulic tortuosity,
y
d In the low Reynolds number regime, flow through a
hλi
- porous matrix is governed by Darcy’s law that links the T = , (3)
u L
fluid flux (discharge per unit area) q with the applied
l
f pressure gradient ∇P by the linear relation
where hλi is the average length of the fluid paths and
.
s L is the geometrical length of the sample [1]. Using the
c k
i q =− ∇P, (1) tortuosity, Kozeny’s relation (3) can be generalized to
s µ
y
φ3
ph where µ is the dynamic viscosity of the fluid and k is a k =c0T2S2. (4)
[ proportionality constant known as permeability [1]. To
a large extent, the proper description of the fluid flow Byfittingexperimentaldata,CarmanconcludedthatT2
1
through a porous medium depends on precise relations is a constant factor (≈5/2) over a wide range of porosi-
v
between the physical properties involved such as perme- ties. Later it was found that T2 does vary with φ, and
6
ability and porosity (φ). In particular, much attention can be as large as 50 for low porosity media [6, 7].
1
3 has been paid to deriving relations between k and φ [2]. Furthermore,itwasrealizedthatelongationofstream-
3 In 1927Kozeny developeda simple capillarymodel for a linesnotonlyaffectsthehydraulicdischarge,butalsome-
. porous medium, and proposed the relation diates other types oftransportphenomena inthe porous
1
0 medium. Thisresultedinintroducingseveraldistinctive,
8 φ3 experimentally measurable tortuosities obtained from a
k =c , (2)
0 0S2 particular transport process, leading to diffusive [8, 9],
: electrical [9, 10, 11] and acoustic [11] tortuosity defini-
v
i whereS isthespecificsurfaceareaandc0 isadimension- tions. There were also further theoretical attempts to
X
less Kozenyconstantthat depends onthe channelgeom- define tortuosity [1, 5, 12]. However, all these tortuosi-
r etry[1]. Unfortunately,Kozeny’sformulaisnotuniversal ties, in general, differ from each other. Except for some
a
anddoes notholdfor complicatedporousgeometries[1]. very simple models [5, 6, 13], there is no clear consen-
For example, it does not take into account pore connec- sus on the relation between these definitions. Among all
tivity and the fact that the specific surface area can be thesedefinitions, the oneexpressedinEq.(3)isnotonly
increasedto an arbitrarily large value by removing some thesimplest,butalsowidelyadoptedintheoreticalstud-
ofthematerialtoroughentheporousmatrixsurfaceina ies, for it ties tortuosity with the underlying geometry
fractal-likemanner. Onapurelyphysicalgroundnamely, and topology of the porous medium.
one would expect that removal of the material from a It has been known since long that flow through a
porous matrix would increase its permeability, whereas porous medium depends on many factors such as poros-
Kozeny’s formula predicts just the opposite [3]. ity, tortuosity, granule shape and size distribution, satu-
One of the most widely accepted attempts to general- ration,Reynolds number, etc. For proper understanding
ize relation (2) was proposed by Carman [1, 4, 5], who of transport phenomena in porous media, however, it is
noticed that the streamlines in a porous medium are far essential to depart from simple systems with a limited
frombeingcompletelystraightandparalleltoeachother. number of well-defined control parameters. Therefore,
2
in this paper we investigate the hydraulic tortuosity (as
defined in Eq. (3)) in a creeping flow through a porous
region constructed by a two-dimensional lattice system
withauniform,randomlydistributedandfreelyoverlap-
ping solid squares. This model, firstused by Koponenet
al. [12], is simple enough to allow a numerical solution
with the advantage of having porosity (φ) as the only
control parameter. It is particularly suitable for study-
ing tortuosity-porosityrelation, as at high porosities the
streamlines are almost straight, whereas for low porosi-
ties they become wiggly. a) φ=0.5 b) φ=0.8
Using the Lattice Gas Automata (LGA) method, Ko-
FIG. 1: An example of two 800×800 porous matrices con-
ponen et al. [12] solved the flow equations for a porosity
structedbyrandomlyplacedandfreelyoverlappingrectangles
range of φ∈[0.5,1], and concluded that
of size 10×10 for two different porosities φ.
T =p(1−φ)+1, (5)
wherepisafittingparameter. However,latertheyfound aligned with the y-axis of an x-y Cartesian coordinate
therelationnotbeingconsistentwiththeresultsobtained system to model the gravity. Following previous works
for the porosity range φ ∈ [0.4,0.5], and suggested [14] [12,14],periodicboundaryconditionshavebeenimposed
to replace relation (5) by in both directions.
Two examples of such porous systems are depicted in
(1−φ) Fig. 1. The dark areas represent fixed solid obstacles,
T =p +1, (6)
(φ−φ )m while the white part is occupied by the fluid.
c
in which φ ≈ 0.33 is the percolation threshold while p
c
andmaresomeempiricalparameters. Still,asanad hoc
B. Numerical techniques
formulawithtwofittingparameters,thisrelationcannot
be consideredauniversallaw. Inaddition,the dataused
Numerical solution of the model defined above con-
to derive(6)suffer fromsystematicerrors,asneither the
sistsoffivemainsteps: (i)generationofaporousmatrix
impact of a finite system size nor the effect of the space
of a known porosity; (ii) solving the flow equations in
discretization were taken into account.
the low Reynolds number regime; (iii) finding the flow
The aimofthis paper is to caryouta numericalsimu-
streamlines; (iv) determining the tortuosity of the flow;
lationforanalysingthetortuosity—porosityrelationina
(v) error analysis.
systemoffreely overlappingrectangleswitha highaccu-
racy. In addition, we provide a simplified algorithm for
T calculationwithoutthe needforimplementing compli-
1. Construction of the porous matrix
cated, weighted averages of streamline lengths.
The structure of the paper is as follows. Section II
specifies the model and the numerical techniques used. Aporousmatrixofagivenporosityφcanbegenerated
Special attention is paid to the description of the non- using the method of [12, 14]. Starting from an empty
trivial numerical technique for the tortuosity. Main re- system,solidsquaresareaddedatrandompositionsuntil
sults, including a detailed numerical error and finite-size the desired porosity has been reached. The porosity is
analysis are provided in Sec. III. Finally, the results are calculated as the fraction of unoccupied lattice nodes.
discussed in Sec. IV.
2. Lattice Boltzmann method for solving flow equations
II. MODEL
To solve the flow equations, we applied the Lattice
A. General description
Boltzmann Model (LBM) [15] with a single relaxation
time collision operator [16]. This method proved useful
The system of interest consists of a square lattice in microscopic model simulations of flow throughporous
(L×L) in which a number of solid squares (a×a lat- mediaforvariousconditionsandflowregimes[2,17]. Itis
tice units)havebeenplacedatrandomlocationsto form anumericaltechniquethatrestsontheBoltzmanntrans-
a porous matrix (1 ≤ a ≪ L). The squares are fixed port equation discretized both in time and space, and is
in space but free to overlap. The only restriction is that expressed in terms of the velocity distribution functions
theirsidesmustcoincidewiththeunderlyinglattice. The n in the form
i
remaining, void space is filled with a liquid. The con-
stant, external force imposed on the porous medium is nt+1(r+c )=nt(r)+∆t(r)+F , (7)
i i i i i
3
A(cid:13) A(cid:13) A(cid:13) A(cid:13)
A(cid:13) A(cid:13)
A(cid:13) A(cid:13) A(cid:13) A(cid:13)
A(cid:13) A(cid:13) C(cid:13) C(cid:13) B(cid:13) B(cid:13)
A(cid:13) C(cid:13) B(cid:13)
A(cid:13) A(cid:13) C(cid:13) C(cid:13) B(cid:13) B(cid:13)
B(cid:13) B(cid:13) B(cid:13) B(cid:13) D(cid:13) D(cid:13) D(cid:13) D(cid:13)
B(cid:13) B(cid:13) D(cid:13) D(cid:13)
B(cid:13) B(cid:13) B(cid:13) B(cid:13) D(cid:13) D(cid:13) D(cid:13) D(cid:13)
D(cid:13) D(cid:13) D(cid:13) D(cid:13)
D(cid:13) D(cid:13)
D(cid:13) D(cid:13) D(cid:13) D(cid:13)
a)(cid:13) b)(cid:13)
FIG. 2: A schematic view of the system before (a) and after
(b) grid refinement. In this examplethe system size is 4×6,
the porous matrix is made of three quads of size 2×2 (A, B
and D), and therefinement level kref =2.
where i = 0,...,8 identifies lattice vectors c , t is an
i
integer time step, r denotes a lattice node, ∆t(r) is the
i
collisionoperatoratr,andF representsthei-thcompo-
i
nent of the external force. We used the time unit equal
to the relaxation time, which yields the kinematic vis-
cosity ν = 1/6 [15]. This, in turn, simplifies the form
of the collision operator ∆t(r) = neq(r) − nt(r) with
i i i
neiq(r) being the equilibrium value of ni at r. The ex- FIG. 3: Velocity magnitudes squared (u2 = u2x+vy2) calcu-
ternal force was taken into account using a method of lated on a600×600 numerical grids, which corresponds to a
Ref. [18]: half ofthe momentum was transferreddirectly 200×200lattice(refinementlevelkref =3). Thesquareblock
intothe equilibriumdistributionfunctionduringthe col- sizeswere10×10(i.e.,30×30afterrefinement)andtheporos-
lision step, whereas the other half was included into the itywasφ=0.65. Periodicboundaryconditionswereassumed
transport equation. Because we were interested in the in both directions. The grey squares in the figure represent
solution of a slow, laminar flow, we utilized the equilib- thesolid part ofthemedium,whereastheblackregion isthe
rium distribution function neq linearized in the velocity porespace open to fluid flow.
i
as
neq =w ρ[1+3(u·c )], (8) no-slipboundaries,secondorderaccuratesolutions,both
i i i
in space and time, could be achieved [15]. An example
in which u is the macroscopicvelocity vector and w are
i of the velocity field calculated with this method for a
some weighting coefficients that depend on the lattice
low-porosity medium is shown in Fig. 3.
structure and dimension [19, 20].
One problem with the LBM method is that it is inca-
pable of resolving the macroscopic Navier-Stokes equa-
3. Flow streamlines
tions for channels narrower than about 4 lattice units
[15]. This limitation becomes particularly important at
low porosities,for whichthe number ofvery narrowpas- After obtaining the velocity field u at eachgrid point,
sages increase enormously. To bypass this problem, a streamlines could be obtained by solving the equationof
standardnumericalmeshrefinementprocedurewasused. motionforthe trajectoriesr(t)ofmasslessparticles[21],
Starting from the original lattice taken to generate the
porousmatrix,eachofitsL2 elementaryquadsweresub- dr =u(r). (9)
divided into k ×k smaller quads with k =1,2,... dt
ref ref ref
being the refinement level. The resulting k L×k L
ref ref The u(r) values of points lying between two grid nodes
computational grid of vectors r in Eq. (7) will then be
were obtained using bilinear interpolation. Due to com-
formed from the centers of the small quads. With this
plex boundary conditions and extreme velocity differ-
choice, the identification of the interface between the
ences onthe grid,the 4-thorderRunge-Kuttaalgorithm
porous matrix and the free space is facilitated.
with adaptive time stepping was used [22].
A schematic picture of the refinement procedure is
shown in Fig. 2. Note that the refinement effectively in-
creases the number of the lattice nodes between any two
4. Tortousity
points by the factor k , and that the smallest channel
ref
width is k +1 lattice units.
ref
After initialization, the LBM computational loop of The tortuosity is defined by Eq.(3) as the ratio of the
advection and collision continued for 5×103 time steps. averagelengthofallparticlepathlinespassingthrougha
By using the mid-grid bounce-back rule applied to the givencross-sectionduringaunittimeperiodtothewidth
4
of the sample [1] leading to
T = 1 RAuy(x)λ(x)dx, (10)
L u (x)dx
RA y
in which A is an arbitrary cross-section of the system
paralleltothexaxis,λ(x)isthelengthofthestreamline
cutting Aatx, andu (x)is the y componentofthe trial
y
particle velocity at x.
The integrals in (10) have been obtained in the liter-
ature either by the Monte Carlo integration [12, 14, 23]
or by direct quadratures [6]. In the former method, the
lengthsofthestreamlinespassingthroughrandomlycho-
sen points within the pore volume are averaged using
proper weights. In the latter method T is approximated
by the relation
1 Pjuy(xj)λ(xj)∆xj
T ≈ , (11)
L u (x )∆x
Pj y j j
where ∆x = x − x are discretization intervals of
j j+1 j
A. In principle, both approaches should yield the same
FIG. 4: Streamlines generated with the constant-flux con-
results, but both can be easily misused. For exam- straint (12) for the same system as in Fig. 3 (N = 59). The
ple, some researchers used the Monte Carlo integration horizontal line represents thecross-section y=L/2 on which
withstreamlinespassingthroughpointschosenrandomly the initial points xj were chosen. Only those streamlines are
from a uniform distribution over the whole pore space shownforwhichλ(xj)couldbedeterminednumerically. The
[12, 23], some others calculated streamlines cutting all arrows point at regions discussed in thetext.
lattice nodes [14], whereas others recorded all stream-
linescrossingeverylatticenodealongachoseninletplane
[6, 13]. However, such ‘uniform’ approaches are not co-
herent with the reality of low porosity systems, in which are two types of such regions: dead-end pores and vol-
transport is mostly carried out only through few ‘con- umes where the fluid stream is split by an obstacle (or
ducting’ channels (cf. Fig. 3). Consequently, the sums merges behind it). An example of xj located in a dead-
in (11) contain most probably many terms of practically end pore is shown by arrow (a) in Fig. 4. Arrow (b)
negligible magnitudes. To avoid this problem, we used points at xj corresponding to an incomplete streamline
Eq. (11) with a constant-flux constraint between two that flows from a region (c) where the stream merges
neighboring streamlines, after being split by an obstacle.
To bypass this problem,in calculationof the sum (13)
u (x )∆x =const. (12)
y j j
only those streamlines were taken into account, whose
lengths λ canbe determined. The errorinduced by this
With this choice, Eq. (11) immediately simplifies to j
procedure is discussed in Sec. III.
N
1 1
T ≈ Xλ(xj), (13)
LN
j=1
where N is the number of the streamlines generated.
Note that all terms in this sum are of the same order 5. Error analysis
of magnitude.
Thus, to calculate T, a horizontal cross-section A is
chosen. Next, the coordinates of the initial points x are Tortuosityvalues(T)calculateddirectlyfrom(13)con-
j
determinedusing(12),andthecorrespondingstreamlines tain errors arising from different sources. While statisti-
are found by solving (9) in both directions until the so- cal errors result from randomness in the porous matrix,
lutions hit the system edges (y = 0 and y = L). Finally discretization errors appear when approximating the in-
their lengths are plugged into (13). tegralsin(10)byfinitesums,andwhensolvingflowequa-
It should be noted that not all streamlines passing tions by discrete lattices. Finite-size errorscouldemerge
through x cut both horizontal edges y = 0 and y = L also as a consequence of approximating a macroscopic
j
(see Fig. 4). This may happen if a streamline passes system with a microscopic model. Details of the error
througharegionwithextremelylowfluidvelocity. There analysis are addressed in the next section.
5
1.6
1.5
1.4
T 1.3
1.2
1.1
f = 0.5
f = 0.7
1
0 50 100 150 200 250 300 350 400 450
L
FIG. 6: Tortuosity T as a function of the system size L for
a=10 and φ=0.7 and 0.5, averaged over M =30 samples.
The solid line is the best fit calculated with Eq. (14). The
error bars represent thestandard error of themean.
continuous (Fig. 5c). Moreover, the problem of finding
the coordinates of discontinuity points is numerically ill-
conditioned. Thesetwofactorsgreatlycomplicatethede-
terminationoftheenumeratorinEq.(10),andintroduce
an additional source of errors in (13). For x near a dis-
j
FIG.5: Thebasicquantitiesthataffecttortuositycalculation
continuity point, even small numerical errors may result
byEq. (10): (a) velocity component uy(x);(b) local tortuos- inasignificantjumpinλ(x ). Twocountermeasureswere
ity τ(x) = λ(x)/L; (c) the product uy(x)τ(x); (d) the ratio takento reducethe impactjofthis phenomenon,whichis
of the minimum to maximum speeds along streamlines. All
closely related to the problem of “missing streamlines”
quantities were determined for the system shown in Fig. 4
withthecross-sectiony=L/2andtheuniformdiscretization discussed in Sec. IIB. First, a check is made to find
∆xj = L/N (N = 1200). The tortuosity for this system is out how T calculated from (13) varies with N. Here, an
T ≈1.45. optimal value of N ≈L was found. Second, the tortuos-
ity was always calculated as an average over 8 different
cross-sections. This approachnot only reduced the error
III. RESULTS resulting fromapproximating(10) by (13), but alsogave
someestimationonitsmagnitude. Theerrorswerefound
tobemaximumforlowporosities,butevenforφ=0.45,
To begin te discussion, the structure of the integrands
the relative error was less than 0.5%.
in (10) is examined. Figure 5 shows u (x), the local tor-
y The large number of discontinuities in τ(x) implies
tuosity τ(x)=λ(x)/L, their product u (x)τ(x), and the
y that the fluid velocity along a typical streamline may
ratio of the minimum to maximum trial particle speeds
vary by many orders of magnitude. This is shown in
for the same system as in Figs. 3 and 4. The data for
Fig. 5d, in which the ratio of the minimum to maximum
this figure were obtained from the streamlines originat-
fluid speed (u /u ) alongthe streamline cutting the
ing at the cross-section y = L/2. As expected, the ve- min max
cross-section y = L/2 at x is plotted. In this particular
locity profile is continuous and piecewise-differentiable,
case u /u ≤ 0.26 and drops to 0 at all positions
and partially resembling that of a Poiseuille (parabolic) min max
whereτ(x) became discontinuous. Acomparisonofpan-
flow. The negative value of u near x = 160 indicates
y els(a)and(d)inFigure5showsthatstreamlinespassing
thatsomestreamlinesarecuttingtheinitialcross-section
through a region with relatively high fluid velocity will
many times. This effect must be taken into account to
likely hit regions where the fluid is almost stagnant. For
avoidmultiple countsofthe samestreamlineinEq.(13).
thisreasonitisessentialthatEq.(9)besolvedwithanu-
In contrast to u , the local tortuosity τ is a discontin-
y mericalmethod that uses variable step lengths and local
uous function of x (Fig. 5b). Eachjump in the τ(x) plot
error control.
corresponds to the stream flow splitting into (or merg-
Next, finite-size effects were analysed. Figure 6 shows
ing from) two parts upon meeting an obstacle. Thus,
the dependency of the tortuosity T on the system size
for a finite-size system, τ(x) is a piece-wise continuous
L for two porosities φ = 0.7 and φ = 0.5 averaged over
function with a certain number of discontinuities, equal
M =30 samples. The lines represent fits to
to the number of ,,islands” existing in the porous do-
main. Consequently, the product uy(x)τ(x) is also dis- T(L)=T∞−bexp(−cL), (14)
6
1.6 1.7
1.55 1.6
1.5 1.5
T 1.45 T 1.4
1.4 1.3
1.35 1.2
k =4
ref
1.3 kref=3 1.1
k =1
ref
1
0 50 100 150 200 250 300 350 400 450 0.4 0.5 0.6 0.7 0.8 0.9 1
L f
FIG. 7: The dependency of T on the system size L for φ = FIG. 8: The dependency of the tortuosity T on the porosity
0.45 and three refinement levels kref (symbols). The lines φ. Our data obtained with the LBM method and Eq. (13)
represent best fits tothe function given by Eq.(14). (circles); relation (6) derived for the same system by Koppo-
nen et al. (dashed line); thebest fit toEq. (15b) (solid line).
where T∞, b and c are free parameters. These, as well
as many other fits, not shown here, proved to be suit-
refinements. Hence, the discrepancy can be explained
able enough to estimate the tortuosity of an infinite sys-
as a consequence of finite-size effects and discretization
tem T∞. This procedure also enabled us to estimate a
∗ error analysis, which are missing in their study.
characteristic system length L above which T does not
∗ We fited our data to four tortuosity-porosity relations
changesignificantly with L. It turned out thatL ≈200
∗ proposed by other researchers:
for φ ≥ 0.55 and L ≈ 300 for 0.45 ≤ φ < 0.55. In all
casesanalysed,T wasfoundto be anincreasingfunction T(φ) = φ−p, (15a)
of L. Thus, it becomes clear that ignoring finite-size ef-
∗ T(φ) = 1−plnφ, (15b)
fectsandusingL<L wouldleadtoanunderestimation
of T. T(φ) = 1+p(1−φ), (15c)
Following this, the sensitivity of tortuosity measure- T(φ) = [1+p(1−φ)]2, (15d)
ments to the numerical mesh refinement was examined.
Figure 7 depicts the values of T(L) for three refinement where p is a parameter. The first of them was proposed
levels kref = 1,3,4. We found that T∞ significantly de- forthe electrictortuositybyArchie(1942)[24]. Thesec-
pendsonkref onlyforkref ≤3. Thisthresholdvalueisin ond equation (with p = 1/2) was found in a theoretical
accord with the criteria, mentioned above, for the LBM study on diffusivity of a model porous system composed
method to reconstruct the Navier-Stokes equations [15]. of freely overlapping spheres by Weissberg (1966) [25].
Notethatignoringnumericalmeshrefinementwouldlead The same relation (with p ≈ 0.86 and p ≈ 1.66) was
to a significant underestimation of T. This shows that also reported in measurements of the hydraulic tortuos-
narrow passages are a relevant factor that affects trans- ityforfixedbedsofparallelepipedalparticleswithdiffer-
port properties of a porous medium. Numerical mesh entthickness-to-sideratiosbyComitiandRenaud(1989)
refinement is thus particularly important at low porosi- [26]andinrecentmeasurementsofelectricaltortuosityin
ties, where narrow passages are common. fixed beds and suspensions of glass spheres by Barrande
Afterfindingtheminimalrequirementsonthemeshre- et al (2007)[27]. Equation (15c) is an empirical relation
finementlevelk andthesystemsizeL∗,thetortuosity- found for sandy (p = 2) or clay-silt (p = 3) sediments
ref
porosity relation could be determined. For a given by Iversen and Jørgensen (1993) [28]. Finally, equation
φ = 0.45,0.5,...,0.95, a system size of L = L∗ and a (15d), with p = 32/9π ≈ 1.1, was recently obtained in
refinement level of k = 3 were chosen. At each φ, T a model of the diffusive tortuosity in marine muds by
ref
was calculated for M porous matrices, with M ranging Boudreau and Meysman (2006) [29].
from 30 (for φ = 0.95) to 140 (for φ = 0.45), and the Even though we treated p in all these formulas as an
results are shown as circle symbols in Fig. 8. For com- adjustable parameter,only Eq. (15b) gave a satisfactory
parison, we also plotted the best-fit curves obtained for fit for p = 0.80±0.01 (the reduced chi-square statistics
exactly the same system by Koponen et al. [14], see Eq. ≈ 0.8), a value comparable with those of Comiti and
(6). Obviously, the results of Koponen et al. lie signif- Renaud [26] (p ≈ 0.86 and p ≈ 1.66). This fit is plotted
icantly below those obtained by us. This is due to the in Fig. 8 as a solid line.
fact that in the work of Koponen et al., a rather small The fact that our systemobeys Eq.(15b) has a rather
∗
system (L < L ) is considered without numerical mesh interesting and unexpected consequence. As shown pre-
7
viously [14], the specific surface area S in the model of fluidflow,trajectoriesofindividualmoleculesarenotlim-
freely overlapping squares of side a satisfies the relation ited to single “theoretical” streamlines, but are affected
by diffusion at low velocities or turbulence at higher ve-
2
locities.
S =− φlnφ, (16)
R It is also interesting to mention that the problem of
finding T is easiest at high porosities, where nearly each
with R = a/2 denoting hydraulic radius of obstacles.
obstacle constitutes a separate island. Although the
With this, equation (15b) simplifies to
number of discontinuities is very large, in general they
S tend to average out. At low porosities, however, se-
T −1∝R . (17)
vereproblemsmayariseasdiscontinuitiesaremuchfewer
φ
in number (which means no “averaging out”) but much
larger in magnitude (which results from increasedisland
IV. DISCUSSION AND CONCLUSIONS sizes). With this study we have demonstrated how sen-
sitive tortuosity computations are to finite-size effects,
discretization errors and large variation of fluid speed
As shown by the present study, obtaining hydraulic
alongstreamlines. Thesystemsizemustbelargeenough
tortuosityfromnumericalsimulationscontainsmanyhid-
to ensure development of chaotic “splitting and merg-
den problems, which may lead to incorrect conclusions.
ing” flows, that are characteristic of real granular sys-
When a fluid stream hits an obstacle, it splits and then
tems. Discretization errors creep into the system from
merges, causing a discontinuity in streamlines. The
severalplaces,mostnotablyinnarrowchannels,andcan
bounding streamline of each obstacle separates the two
be avoided by numerical mesh refinement. Our results
splitting (or merging) streams. The location of such
concerning large fluid velocity variations along stream-
streamlines is a priori not known, and the problem of
lines are a clear indication for revising those tortuosity
finding the streamlines within those regions is numeri-
definitions which assume a constant fluid velocity along
cally ill-conditioned. If the system is sufficiently large,it
a streamline [6, 13]. They also show that numerical de-
is inevitable that majority of streamlines pass through
termination of streamlines requires using advanced nu-
such “ill-conditioned” regions. Moreover, the velocity
merical integrators with adaptive step lengths and local
magnitudealongthestreamlinescanvarybymanyorders
error control.
of magnitude.
These problems resemble those encountered in an- When streamlines are generated using the constant-
othertypeofcomputersimulations,moleculardynamics, flux constraint, the tortuosity can be calculated simply
wherediscontinuities(aswellaslargevariationsinveloc- as an average over the streamline lengths. This method
ity) are caused by collisions. After several such events, reduces the computation errors and does away the need
the computer-generated particle trajectories have liter- for determining the streamlines in dead-end pores.
ally nothing to do with the exact solutions. Still, molec- The numerical data presented in this study were
ular dynamics is one of the most successful methods of found to be in good agreement with those generated by
computer physics. The reason is that physical quanti- Eq. (15b), obtained from previous experimental studies.
ties never depend on exact trajectories of the individual However,thisrelationcannotbeappliedclosetotheper-
particles—itissufficienttoensurethatthesolutionkeeps colation threshold, where tortuosity diverges. Note also
ona constantenergysurface. This analogymakesus be- thatEq. (15b)wasfoundtodescribebothhydraulicand
lieve that despite all difficulties, hydraulic tortuosity is diffusivetortuosities,i.e.quantitiesthatarecertainlycor-
a well-defined quantity that can be reliably obtained by related, but in a way that has not been well established
numerical methods. All simulations we performed with yet. It is not clear whether our findings reflect equiva-
different cross-sections,different numbers ofstreamlines, lence of these quantities, or whether they can be consid-
different choices of streamline starting points, different ered as a coincidence.
numericalordinary differential equation solvers,resulted We also found that in the model of freely overlapping
inalmostthesamenumericalvalues. Thisimpliesthat— squares,a verysimple relation(17) holds between tortu-
just as in molecular dynamics—small local errors,which osity, porosity and the specific surface area. This equa-
areunavoidableincomputersimulations,areofmarginal tionis closelyrelatedto Eq.(15b), andit expectedto be
importance. This corresponds to the fact that in real valid for all systems, for which Eq. (15b) holds.
[1] J. Bear, Dynamics of fluids in porous media (Elsevier, [3] A. W. J. Heijs and C. P. Lowe, Phys. Rev. E 51, 4346
NewYork, 1972). (1995).
[2] A.Koponen,D.Kandhai,E.Hell´en,M.Alava,A.Hoek- [4] P. C. Carman, Trans. Inst.Chem. Eng. 15, 150 (1937).
stra,M. Kataja, K.Niskanen,P.Sloot, andJ.Timonen, [5] M.B.Clennell,GeologicalSociety,London,SpecialPub-
Phys.Rev.Lett. 80, 716 (1998). lications 122, 299 (1997).
8
[6] M. A.Knackstedtand X.Zhang, Phys.Rev.E 50, 2134 [19] R. Verberg and A. J. C. Ladd, Phys. Rev. E 60, 3366
(1994). (1999).
[7] B. P. Boudreau, Geochimica et Cosmochimica Acta 60, [20] X. Heand L.-S.Luo, Phys. Rev.E 56, 6811 (1997).
3139 (1996). [21] F. H. Harlow and J. E. Welch, Phys. of Fluids 8, 2182
[8] Y. Nakashima and T. Yamaguchi, Bulletin of the Geo- (1965).
logical Surveyof Japan 55, 93 (2004). [22] W.H.Press,B.P.Flannery,S.Teukolsky,andW.T.Vet-
[9] A.A. Garrouch, L. Ali, and F. Qasem, Ind.Eng. Chem. terling, Numerical Recipes: The Art of Scientific Com-
Res. 40, 4363 (2001). puting (Cambridge University Press, Cambridge (UK)
[10] P.B. Lorenz, Nature189, 386 (1961). and New York,1986).
[11] D. L. Johnson, T. J. Plona, C. Scala, F. Pasierb, and [23] P.Alam,T.Byholm,andM.Toivakka,NordicPulpand
H.Kojima, Phys.Rev.Lett. 49, 1840 (1982). Paper Research Journal 21, 670 (2006).
[12] A. Koponen, M. Kataja, and J. Timonen, Phys. Rev. E [24] G. Archie, Transactions of the American Institute Min-
54, 406 (1996). ing, Metallurgical and Petroleum Engineers 146, 54
[13] X.ZhangandM.A.Knackstedt,Geophys.Res.Lett.22, (1942).
2333 (1995). [25] H. L. Weissberg, J. Appl.Phys. 34, 2636 (1963).
[14] A. Koponen, M. Kataja, and J. Timonen, Phys. Rev. E [26] J. Comiti and M. Renaud, Chem. Eng. Sci. 44, 1539
56, 3319 (1997). (1989).
[15] S. Succi, The Lattice Boltzmann Equation for Fluid Dy- [27] M. Barrande, R.Bouchet, and R. Denoyel, Anal. Chem.
namicsand Beyond (ClarendonPress,NewYork,2001). 79, 9115 (2007).
[16] P.L.Bhatnagar, E.P.Gross, andM.Krook, Phys.Rev. [28] N. Iversen and B. B. Jørgensen, Geochimica et Cos-
94, 511 (1954). mochimica 57, 571 (1993).
[17] C. Pan, M. Hilpert, and C. T. Miller, Phys. Rev. E 64, [29] B. P. Boudreau and F. J. Meysman, Geology 34, 693
066702 (2001). (2006).
[18] Z. Guo, C. Zheng, and B. Shi, Phys. Rev. E 65, 046308
(2002).
