VOLUME 80, NUMBER 4 PHYSICAL REVIEW LETTERS 26 JANUARY 1998

Permeability of Three-Dimensional Random Fiber Webs

A. Koponen,1 D. Kandhai,2 E. Hellén,3 M. Alava,3,4 A. Hoekstra,2 M. Kataja,1 K. Niskanen,5

P. Sloot,2 and J. Timonen1
1
Department of Physics, University of Jyväskylä, P.O. Box 35, FIN-40351 Jyväskylä, Finland
2
Department of Mathematics, Computer Science, Physics and Astrophysics, University of Amsterdam,
Kruislaan 403, NL-1098 SJ Amsterdam, The Netherlands
3
Laboratory of Physics, Helsinki University of Technology, P.O. Box 1100, FIN-02015, HUT, Finland
4
Nordita, Blegdamsvej 17, DK-2100 Copenhagen, Denmark
5
KCL Paper Science Center, P.O. Box 70, FIN-02151 Espoo, Finland
(Received 25 August 1997)
We report the results of essentially ab initio simulations of creeping flow through large three-
dimensional random fiber webs that closely resemble fibrous sheets such as paper and nonwoven
fabrics. The computational scheme used in this Letter is that of the lattice-Boltzmann method and
contains no free parameters concerning the properties of the porous medium or the dynamics of the
flow. The computed permeability of the web is found to be in good agreement with experimental
data, and confirms that permeability depends exponentially on porosity over a large range of porosity.
[S0031-9007(97)05087-4]

PACS numbers: 47.55.Mh, 47.15.Gf

Fluid flow in porous media plays an important role in mented by the recent introduction of new flow-simulation
a wide variety of technical and environmental processes algorithms that are particularly suited for parallel com-
such as oil recovery, paper manufacturing, and spread of puting. Among these are the lattice-gas-automaton [7,8]
hazardous wastes in soils. The single-phase creeping flow and lattice-Boltzmann method [9–11], which have al-
through a porous substance in a gravitational field g is ready been applied to a wide class of fluid-flow problems.
well described by Darcy’s law [1,2] which, in this case, These methods are particularly useful in complex and ir-
can be written in a form regular geometries [12 –16].
k Despite the numerous experimental and theoretical stud-
q­ g, (1) ies (see, e.g., [5] for a comprehensive review), permeabil-
n
where q is the fluid flux through the medium, n is the ity characteristics of disordered fibrous porous media are
kinematic viscosity of the fluid, and k is the permeability still poorly understood. The existing numerical studies
coefficient that is a measure of the fluid conductivity include those on fluid flow through random arrays of par-
through the porous medium. allel cylinders, suspension of prolate spheroids, and three-
Determination of coefficient k for each particular dimensional regular fiber networks, which all neglect the
substance has been a long-standing problem. The ex- disorder typical of real 3D fiber webs [17–19].
perimental methods that have been used to this end vary In this Letter we report the results of three-dimensional
from rather straightforward measurements [3–5] to more lattice-Boltzmann simulations of creeping flow through
sophisticated approaches, which utilize, e.g., mercury large random fiber webs, and compare these results with
porosimetry, electrical conductivity, nuclear magnetic previous experimental, analytical, and numerical results
resonance, or acoustic properties of the medium [6]. for various fibrous materials.
Theoretical methods typically rely on analytical models The model web structures were constructed in dis-
based on simplified pore geometries, which allow solution cretized space using a recently introduced growth algo-
of the microscopic flow patterns [1], or on more advanced rithm [20]. Within this algorithm, fiber webs are grown
methods that statistically take into account the structural by sequential random deposition of flexible fibers of rect-
complexity of the medium [1,2]. angular cross section on top of a flat substrate (the xy
Numerical simulations are often used to connect the- plane). Each fiber was randomly oriented either in the x
ory with experiments. Realistic three-dimensional flow or y direction, and was then let to fall in the negative z
simulations in complex geometries are, however, very de- direction until it made its first contact with the underlying
manding in terms of computing power. Until recently structure. After this it was bent downwards without de-
this approach has thus been hampered by the necessity structing the structure, and subject to the constraint
of major simplifications in the pore structure or flow jzi 2 zj j # F . (2)
dynamics. New techniques based on massively parallel
computers and increased single processor capabilities Here zi and zj are the elevations of the fiber surface above
have now made 3D ab initio simulations of realistic two nearest-neighbor cells i and j, and F is an effective
flow problems feasible. This development is further aug- fiber flexibility. Notice that for long fibers, the resulting

716 0031-9007y98y80(4)y716(4)$15.00 © 1998 The American Physical Society

VOLUME 80, NUMBER 4 PHYSICAL REVIEW LETTERS 26 JANUARY 1998

web structure, porosity, and contact area of fibers depend grid resolution. Resolution sensitivity is caused by the
only on F [21]. A large F produces a dense web while a bounce-back boundary condition applied to the solid-fluid
small F leads to a more porous structure. interface, and by Knudsen effects in small pores [8,11,14].
In order to construct structures that are homogeneous in These effects are viscosity dependent, and they determine
the z direction, the samples used in the simulations were the minimum size of obstacles and pores that can be used
extracted from inside of thicker webs. Ten grid layers of in simulations.
void space were then added on the top of samples, and In order to find an acceptable grid resolution in different
the system was made periodic in all directions. In Fig. 1 porosity regions, we first made a series of test runs using
we show a sample created by this algorithm. It is evident fibers of aspect ratio 10, i.e., of size wF 3 wF 3 10wF ,
that the produced structures closely resemble those of, e.g., with wF ­ 5, 10, and 20 grid units. In these tests the
paper [22] and nonwoven fabrics (restriction to the x and size of the simulation sample was Lx 3 Ly 3 Lz ­
y directions can be relaxed and does not play an important 20wF 3 20wF 3 10wF . The relaxation parameter t of
role here). the LBGK collision operator was varied from 0.668 to 2
Simulations of fluid flow through the web in the z corresponding to the dynamic viscosity n ­ s2t 2 1dy6
direction were done using the 19-link LBGK (lattice- [9] ranging from 0.056 to 0.5 (in grid units).
Bhatnagar-Gross-Krook) model [9]. Flow was induced In Fig. 3 we show the computed permeability as a
by applying a “gravitational” body force on the fluid [14]. function of viscosity n for two test systems of different
This was accomplished simply by adding at each time step porosities. The nearly linear dependence on n of perme-
a fixed amount of momentum in the negative z direction ability [14] is clearly seen in this figure. For f ­ 0.67
to all “particles” within the pore space. In Fig. 2 we [Fig. 3(a)], the result is already almost independent of the
show the simulated stationary velocity field for a flow in grid resolution for the smallest value of viscosity. For
the z direction through the highly inhomogeneous sample high porosities f (with n ­ 0.056), resolution wf ­ 5
shown in Fig. 1. can thus be used. For a smaller f finite-size effects be-
It is evident that there are large fluctuations in the come more pronounced. Comparing the simulated perme-
velocity field reflecting the variations in the local porosity abilities for f ­ 0.39 [Fig. 3(b)] at n ­ 0.056, we can
of the sample. The average velocity kyz l shown in Fig. 2 conclude that wF ­ 10 is satisfactory for low porosities.
is kyz l ­ 0.001 42 (in grid units), with a standard deviation In the actual permeability simulations fibers of aspect
of Dyz ­ 0.001 28. These fluctuations, which are inherent ratio 20, i.e., of size wF 3 wF 3 20wF , and a sample
in random porous structures, will affect the permeability, of dimensions 80wF 3 80wF 3 10wF were used. (In
except at very high porosities, such that it is expected to Fig. 1 we show an example of the webs for which
become higher than that for regular arrays of pores [5]. simulations were performed.) The flexibility parameter F
This effect will be seen in the results given below. was varied from 0 to 3 which corresponds to porosities f
In a stationary state the body force applied to the fluid ranging from 0.95 to 0.42. In the simulations t ­ 0.668
is completely canceled by viscous friction forces due to (n ­ 0.056) and two different values for wF were used.
fibers. Once the total flux of fluid through the sample For f . 0.6, wF ­ 5 was used, and the size of the
and the viscosity of the fluid are known, the permeability simulation lattice was 400 3 400 3 60 grid points. For
of the sample can then be determined from Darcy’s law f , 0.6, wF ­ 10 and a lattice of 800 3 800 3 110
Eq. (1). It is well known that the permeability of a porous grid points were used. For these discretizations, the
medium determined by the present method depends on

FIG. 2(color). The velocity field of fluid flow through the

FIG. 1(color). A fiber-web sample constructed with the depo- fiber web shown in Fig. 1. Bright colors indicate high fluid
sition model. The porosity of the web is 0.83. velocity.
717
VOLUME 80, NUMBER 4 PHYSICAL REVIEW LETTERS 26 JANUARY 1998

FIG. 4. The calculated dimensionless permeability kya2 as

FIG. 3. Calculated dimensionless permeability kya2 as a a function of porosity (black triangles). Open squares and
function of viscosity n for two test samples with porosities circles show the experimental results for fibrous filters [5] and
(a) f ­ 0.67 and (b) f ­ 0.39. Here a ­ wF y2 is the hy- compressed fiber mats [3,5], respectively. Curve (1) is the
draulic radius of the fibers. The fiber widths (grid resolutions) analytical result for cubic lattice given in Ref. [5], curve (2)
are wF ­ 5 (squares), 10 (circles), and 20 (star). is the numerical result for an fcc lattice from Ref. [19], and
curve (3) shows the result of a fit with the Kozeny-Carman
relation, Eq. (3).
estimated finite-size errors of the simulated permeabilities
were less than 15%.
Simple dimensional analysis suggests that, for a constant so that there are significant hydrodynamic interactions
body force, the saturation time for reaching the steady state between them. On the other hand, the simple capillary-
2
is proportional to Rpore yn, where Rpore is the characteristic tube model by Kozeny and Carman [1] gives kya2 ~
size of pores. For systems with high porosity, saturation f 3 s1 2 fd22 in this limit so that, as expected, the
time thus tends to become very long. This time can, simulated behavior of ksfd is rather in good agreement
however, be essentially reduced by using an iterative with this model when f & 1.
momentum relaxation (IMR) method in which the applied A fit of the form lnskya2 d ­ A 1 Bf to the rest of
body force is changed during simulation depending on the simulated points gives A ­ 28.53, B ­ 10.4, with a
the value of the friction force. A thorough description very high correlation between the simulated points and the
of the IMR method and a detailed benchmark analysis fitted curve. So far there have been no analytical results
of our parallel lattice-Boltzmann code will be published which would have produced this kind of exponential
separately [23,24]. In the IMR method the number of time behavior at intermediate porosities. It will not hold near
steps required for saturation (with 1% accuracy) varied the percolation threshold at which permeability vanishes.
between 1000 to 8000 depending on porosity, and was This critical region is, however, beyond the present
typically less than half of that needed in the constant- computational capabilities.
body-force method. When 32-bit floating point numbers Experimental results [3,5], shown in Fig. 4 as open
were used, the larger simulation lattice required 5.4 Gbyte circles and squares, conform well with the simulated
of core memory. The simulations were therefore carried points. Notice that there is no free parameter in the
out using 64 nodes (300 MHz CPUs with 128 Mbyte of present model as permeability is scaled by the square of
memory) on a Cray T3E system. The required CPU time the hydraulic radius of the fibers, and the relaxation pa-
was typically between 1 and 4 hours. rameter t is used only to fix the necessary grid resolution.
In Fig. 4 we show the simulated permeability of the The level of agreement is therefore astonishingly high. It
random fiber web as a function of its porosity. In this also shows that the model web used here captures the es-
figure solid triangles denote the simulated values: the sential features of the fibrous filters and compressed fiber
related porosities were given by the flexibility parameter mats used in the experiments.
sFd used. It is evident that there are two distinct features Also shown in Fig. 4 are three curves which are re-
in the simulated ksfd curve. First, it seems to diverge sults of previous analytical [5] [curve (1)], numerical
as expected when f ! 1, and, second, k seems to be an [19] [curve (2)], and semiempirical [1,3,4] [curve (3)]
exponential function of f for a rather wide range of f: considerations. Curve (1) is given by kya2 ­ 2 lns1 2
0.42 & f & 0.85. fd 2 0.931 [with Oss1y lns1 2 fddd], an expression ob-
A fit to the last five points with highest porosities of the tained [5] for a cubic lattice model, and curve (2) results
form kya2 ­ const 3 s1 2 fd2m gives m ­ 1.92. This from a numerical solution [19] for the Stokes flow in a
result shows that even for these rather high porosities face-centered-cubic (fcc) array of fibers. Both these
f & 0.96, the system is not in the true asymptotic region curves are below the simulated points, especially for de-
for which m ­ 1 [25]. The reason for this is that the creasing porosity. Notice that the fcc result also fol-
fibers of the web are still close enough to each other lows an exponential law at intermediate porosities. It is
718
VOLUME 80, NUMBER 4 PHYSICAL REVIEW LETTERS 26 JANUARY 1998

30%–40% below the random fiber-web result for these We thank the Center for Scientific Computing in
porosities, but approaches the latter for small porosities Finland (Sami Saarinen and Tomi Salminen in particular)
since the percolation threshold of the fcc lattice [19] is at for providing computational resources, technical support,
a much lower porosity. and expertise on visualization.
Curve (3) is the Kozeny-Carman expression [1]
k ­ f 3 ycS 2 , (3) [1] J. Bear, Dynamics of Fluids in Porous Media (Dover, New
York, 1972).
where S is the specific surface area of the web, and c is a [2] M. Sahimi, Flow and Transport in Porous Media and
constant in capillary-tube models but is known to depend Fractured Rock (VCH, Weinheim, 1995).
on porosity in fibrous materials [4]. An empirical fit to [3] W. L. Ingmanson, B. D. Andrews, and R. C. Johnson,
measured porosities of these materials gives [3,4] c ­ Tappi 42, 840 (1959).
1
3.5f 3 f1 1 57s1 2 fd3 gys1 2 fd 2 , and we have used [4] S. T. Han, Pulp Pap. Mag. Can. 70, T134 (1969).
this expression to get the curve (3) from Eq. (3). The [5] G. W. Jackson and D. F. James, Can. J. Chem. Eng. 64,
specific surface area S was determined from the surface 364 (1986).
[6] See S. Kostek, L. M. Schwartz, and D. L. Johnson, Phys.
area of the (straight) fibers used to construct the web by
Rev. B 45, 186 (1992); A. H. Thompson, S. W. Sinton,
subtracting the area of the interfiber contacts. Because S. L. Huff, A. J. Katz, R. A. Raschke, and G. A. Gist,
of bending of the fibers this expression gives a lower J. Appl. Phys. 65, 3259 (1989); D. L. Johnson, D. L.
bound for S, and curve (3) is expected to overestimate the Hemmick, and H. Kojima, J. Appl. Phys. 76, 104 (1994),
permeability. This is indeed what happens (cf. Fig. 4). and references therein.
Encouraged by the exponential behavior at intermediate [7] U. Frisch, D. d’Humiéres, B. Hasslacher, P. Lallemand,
porosities of ksfd, we have made an interpolation formula Y. Pomeau, and J.-P. Rivet, Complex Syst. 1, 649 (1987).
that connects this behavior with the right asymptotics in [8] D. H. Rothman and S. Zaleski, Lattice Gas Cellular Au-
the limit f ! 1. We find that the expression kya2 ­ tomata (Cambridge University Press, Cambridge, England,
AfeBs12fd 2 1g21 with A ­ 5.55, B ­ 10.1 fits all the 1997).
simulated points very well. So far we have no theoretical [9] Y. H. Qian, D. d’Humiéres, and P. Lallemand, Europhys.
Lett. 17, 479 (1992).
arguments to support this very simple form for the
[10] R. Benzi, S. Succi, and M. Vergassola, Phys. Rep. 222,
permeability. 145 (1992).
In conclusion, we used the lattice-Boltzmann method [11] I. Ginzbourg and D. d’Humiéres, J. Stat. Phys. 84, 927
on a massively parallel computer to solve ab initio the (1996).
permeability of a large random 3D fiber web as a function [12] A. Cancelliere, C. Chang, E. Foti, D. H. Rothman, and
of its porosity in a large porosity range. An iterative S. Succi, Phys. Fluids A 2, 2085 (1990).
momentum relaxation method was used to considerably [13] A. J. C. Ladd, J. Fluid Mech. 271, 285 (1994); A. J. C.
reduce the computing time needed for reaching stationary Ladd, J. Fluid Mech. 271, 311 (1994).
flows. The web samples used were constructed by a [14] B. Ferréol and D. H. Rothman, Transp. Porous Media 20,
growth algorithm [20] that produces random structures 3 (1995).
of flexible fibers, closely resembling those of, e.g., paper [15] N. S. Martys and H. Chen, Phys. Rev. E 53, 743 (1996).
[16] J. Kaandorp, C. Lowe, D. Frenkel, and P. Sloot, Phys.
sheets and nonwoven fabrics. The simulated results were
Rev. Lett. 77, 2328 (1996).
found to be in excellent agreement with experiments [17] C. K. Ghaddar, Phys. Fluids 7, 2563 (1995).
on the materials mentioned above. There were no free [18] I. L. Clayes and J. F. Brady, J. Fluid Mech. 251, 443
parameters in the simulations which would have been (1993).
used for fitting with the experimental results, and also in [19] J. J. L. Higdon and G. D. Ford, J. Fluid Mech. 308, 341
this sense the results reported are ab initio. We found that (1996).
the exponential dependence on porosity of permeability in [20] K. Niskanen and M. Alava, Phys. Rev. Lett. 73, 3475
a wide range of porosities is a generic feature of fibrous (1994).
porous materials, independent of whether they are random [21] E. K. O. Hellén, M. J. Alava, and K. J. Niskanen, J. Appl.
or not. So far there is no theoretical explanation for Phys. 81, 6425 (1997).
this phenomenon, but a simple interpolation expression [22] K. Niskanen, N. Nilsen, E. Hellén, and M. Alava, in The
Fundamentals of Papermaking Materials, edited by C. F.
for ksfd was formulated which seems to give the correct
Baker (Pira Int., Leatherhead, UK, 1997).
behavior in the whole range of porosities above the critical [23] D. Kandhai, A. Koponen, A. Hoekstra, M. Kataja,
region near the percolation threshold. The results reported J. Timonen, and P. Sloot (to be published).
here clearly demonstrate that ab initio simulations of [24] D. Kandhai, A. Koponen, A. Hoekstra, M. Kataja,
flow in complicated structures are possible (see also J. Timonen, and P. Sloot, Comput. Phys. Commun. (to be
Refs. [14,15]), and many open problems in this area can published).
and will now be addressed. [25] G. K. Bachelor, J. Fluid Mech. 52, 245 (1972).

719