- Email: [email protected]
Efficient simulation of flow and transport in porous media
Chemxol EngineeringScience.Vol. 46. No. ‘3, pp. 2225-2233, 1991. Pnnted in Greal Britain
EFFICIENT
am!-2509/91 S3.03+ 0.00 0 1991 Pergamon Press plc
SIMULATION OF POROUS
FLOW AND MEDIA
TRANSPORT
IN
MUHAMMAD SAHIMI + and DIETRICH STAUFFER* Supercomputer Center HLRZ, c/o KFA, D-5170 Jiilich 1, Germany (First receiued
in revised form 20 November 1990)
10 Auyusf 1990; accepted
Abstract-We describe two efficient computer simulation methods for investigating flow and transport processes in porous media. The first method, which used random-walk techniques, is used for estimating effective diffusion coefficients in disordered porous media. The second method, which is based on lattice gases and celIular automata models, is used for calculating the effective permeability of a porous medium. We describe in some detail the computer algorithms that employ these methods. Etticlent use of vectonzation enables us to simulate diffusion in disordered porous media, representedby two- or three-dimensional lattices, with up to IO9sites, and a speed of nearly four steps per microsecond and per Gray-YMP processor. Vectorization of the algorithm based on cellular automata also enables us to simulate flow in a twodimensional model porous medium with more than a million sites. As examples, we study dilksion and flow in stratified and macroscopically heterogeneous porous media and Row through periodic arrays of obstacles. The relation between the effectivediffusivity and permeability of such porous media is-also briefly
discussed.
1. INTRODUCTION The problems of flow and transport through media are of great theoretical and practical
porous
interest (see, e.g., Scheidegger, 1974; Mason and Malinauskas, 1983; and references therein). For example, in the petroleum industry, the reservoir engineers have been dealing with such problems since oil production began. Problems such as drainage and imbibition in soil, multiphase flow through trickle-bed reactors, mercury porosimetry for determining the pore size distribution of a porous catalyst, channeling in packed columns and ground water flow are but a few processes of interest to hydrologists, soil scientists, and chemical and petroleum engineers. Historically, flow and transport processes in porous media have been simulated starting from continuum equations in which the relevant variables are related to conserved quantities, such as local density, local momentum flux, etc. Macroscopic properties, such as the effective transport coefficients, are then defined as averages of the corresponding microscopic or local quantities (see, e.g., Whitaker, 1986; and references therein). The averages are taken over a volume which is small compared to the volume of the system, but large enough for the transport equations to hold when applied to the volume. At every point of a porous medium one uses the smallest such volume and, thereby, generates macroscopic field variables obeying equations such as Fick’s law of diffusion and Darcy’s law of flow in porous media. This approach has been useful for macroscopic modelling and has been used
‘Author to whom all correspondence should be addressed. Present and permanent address: Department of Chemical Engineering, University of Southern California, Los Angeles, CA 90089-1211, U.S.A. z Present address: Institute for Theoretical Physics, Cologne University, D-5CMl Cologne 41, Germany.
extensively in the past. However, this approach is purely phenomenological and provides little insight into how flow and transport processes depend on the morphology of the pore space, i.e. its topology or pore connectivity, and geometry which is usually represented by a pore size distribution. In fact, in order to use such averaged transport equations, transport coefficients must a priori be known. Therefore, if the morphology of the pore space is strongly chaotic or disordered, a continuum approach may not be a quantitative method of studying transport processes in porous media. But when applicable, one usually uses numerical simulations to solve the continuum equations and obtains a few of the quantities of interest. Such numerical methods often involve some kind of discretization of the continuum equations. The second method of studying flow and transport processes in porous media is based on network or lattice models of porous media. One models the medium as a network of bonds, which represent the pore throats of the actual medium, and sites, which represent the pore bodies of the real system. The chaotic morphology of the medium can be incorporated into the network model by assigning randomly-selected effective sizes to the bonds and sites, to mimic the disordered geometry of the medium, and by using networks of random local coordination number, i.e. the number of bonds connected to the same site, to represent the chaotic topology of the pore space. It is often true that, if the coordination number of a regular network and the average coordination number of a random network are the same, the macroscopic properties of the two networks are very close (see, e.g., Jerauld et al., 1984). so that a network of fixed coordination number is often an adequate representation of the topology of the pore space. In a variation of the network models, one starts with a regular lattice and distributes (at random or in a prescribed fashion) 2225
2226
MUHAMMAD SAHIMI and DIETRICH STAUFFER
inclusions in the system in the form of disks, spheres, ellipsoids, or any other desired shape. The effective size of the inclusions can be fixed or a distributed quantity, and flow and transport take place through the channels between the overlapping or nonoverlapping inclusions. This is the model that we use here to study flow in porous media. Sahimi et al. (1990) have recently reviewed such network modelling approaches to the problems of estimating the effective transport properties of porous media. It can also be shown that there is a close analogy between such network or lattice models and the set of discretized equations that one obtains in the numerical simulations of the continuum equations, so that the network models can be justified on this ground also. Recent advances in computer simulation techniques (see, e.g., Binder, 1987, especially Chapter 8) have now made the network or lattice models an efficient tool of investigating flow and transport processes in porous media. In this paper, we describe two highly efficient methods for studying such processes in porous media. The first method is based on randomwalk techniques (see, e.g., Barber and Ninham, 1970) by which one can obtain reliable estimates of diffusivities and conductivities of disordered porous media. The second method is based on cellular automata (CA) which, in the present context, are essentially the discrete solution of the Navier-Stokes equations. We use the CA method for estimating the effective (Darcy) permeability of a porous medium. The random-walk methods of estimating the transport properties are of course well known. However, recent advances in supercomputer technology have made it possible to turn these methods into efficient computational tools, and study flow and transport processes in lattice or network models of porous media which can have a realistic (and often complex) geometry, using very large systems with millions of sites and bonds. The plan of this paper is as follows. In the next section we describe in some detail the random-walk and CA approaches, and we discuss how these methods can be vectorized for use in a vector computer. In section three, we use these methods to study flow and transport processes in stratified porous media, a highly important problem (see, e.g., Matheron and de Marsily, 1980, Bouchaud et al., 1990), and flow through periodic arrays of obstacles, a long-standing problem of considerable interest (see, e.g., Zick and Homsy, 1982). The paper is summarized in Section 4, where we also discuss the application of these methods to more complex problems. 2.
SIMULATION
simulations on vector or parallel computers like a Cray. A primitive example of CA is Prussian “bureaucracy”, where each site at the next time step assumes the state its superior site has at the present time step. For hydrodynamics, such automata are too simple, and actual computer codes for twodimensional systems use 128 or 256 different states for each site, i.e. 7 or 8 computer bits. The CA approximation of the Navier-Stokes equation in two dimensions (Frisch et al., 1986; Doolen, 1991) is based on particles either resting on the sites of a large triangular lattice, or moving with unit velocity on one of the six bonds emanating from each lattice site. If particles hit each other they are scattered according to the laws of momentum conservation. For example, one particle hitting another particle which is at rest, may be scattered such that its direction changes by MY’, and that the previously resting particle moves in a direction inclined at - 60” with respect to the direction of the incoming velocity (Fig. 1). The inverse process is also possible. If a particle hits a solid wall, it is reflected by 180” (no slip boundary condition). The triangular lattice geometry ensures that in this momentum transfer the total momentum vector is conserved. It also approximates the Navier-Stokes equations better than the square lattice in the continuum limit, because of the isotropy of the fourth-order tensor arising in this limiting process. The extension to three dimensions is more complex since no regular three-dimensional lattice is isotropic and, thus, in the continuum limit, one has spurious terms, in addition to those in the Navier-Stokes equation, which are caused by the anisotropy of the lattice. However, there are now several methods of circumventing this difficulty. For example, one can either use disordered neta three-dimensional and topologically work, such as the Voronoi network which is macroscopically isotropic. Alternatively, one can use (d’Humieres et al., 1986) a four-dimensional facecentered hypercubic lattice (FCHC). For this lattice, which has a coordination number of 24, all pairwise symmetric fourth-order tensors are isotropic and, thus, one is in a position to simulate four-dimensional Navier-Stokes equation on such a lattice. One observes that any solution of the four-dimensional equation which does not depend on the fourth dimension is a solution of the three-dimensional equations. This suggests the use of a FCHC lattice which wraps around periodically in the fourth direction. One can actually use a lattice which is only one lattice unit
high-speed
METHODS
Cellular automata (see, Wolfram, 1986, for a collection of reviews) are large lattices where each lattice site can be in one of several discrete states. The state of each site at the next time step is determined completely by the present state of the neighboring sites. Thus both time and space are discrete, and connections are between neighbors only, ideal conditions for
Fig. 1. Example for a transition rule in the cellular automata method. A particle moving with unit velocity hits a particle at rest, and both are scattered such that the total momentum is conserved in the triangular lattice.
Efficient
simulation of Row and transport
long in the fourth dimension and, therefore, have an effectively three-dimensional spatial structure. The disadvantage of this method is that, although the fourth dimension is very thin, the discrete velocities still have components in all directions; therefore, the model is somewhat bit intensive (24 or 25 bits per site). A third approach (d’Humieres et al., 1986) is to use a three-dimensional cubic lattice in which the particles move with speed zero, one and $ (instead of unit velocity in the two-dimensional lattices). This model uses only 19 bits per site. Despite the discrete nature of CA in both two and three dimensions, these models are capable of exhibiting rich macroscopic complexity such as turbulence (d’Humieres and Lallemand, 1986). In the present paper, the porous mednnn is put into a channel with fixed walls. Obstacles are then distributed in between the walls. In the limiting case of unit porosity (with no obstacles) a simple Poiseuille Aow is expected at low enough Reynolds numbers. At the beginning of the simulation, we construct a transition table which tells us how one of the present states (determined by velocity of the incoming particles) is transformed into the next state, determined by the outgoing velocities. This table contains all of the possible states (which, in our simulations, number 2’ = 256). We then start with a biased random selection of molecular directions, where the bias is such that the average local velocities obey the Poiseuille profile, and simulate the scattering of the molecules by the obstacles of which the porous medium consists of, using the transition table we constructed. There can be two-, three- and four-body collisions at a site, which conserve energy and momentum. The rules are such that no new particle is created, and in one time step all particles on the lattice move first to different site, and then undergo collisions at the new site. One iteration consists of updating all lattice sites according to the transition table, and of the order of lo3 iterations are usually needed to reduce the average velocity to a value which no longer diminishes significantly in further iterations. In this paper, data were taken using 3000 such iterations. The CA method can be used for any porosity, but for low values of porosity, one needs a larger number of iterations to achieve stable results. In these simulations, we arrange circular obstacles in regular periodic arrays as shown in Fig. 2. Fluid flow then takes place in between nonoverlapping
Fig. 2. Flow between
a periodic
arrangement
in porous
media
2227
circles, and from the pressure and velocity distributions the effective permeability of the medium is calculated, just as in the case of one single circle or many randomly-placed and partially overlapping circles (Brosa and Stauffer, 1990; Duarte and Brosa, 1990). Various computational implementations of these ideas are discussed in the literature (Doolen, 1991); we follow Brosa and Staufier (1989) who published the innermost FORTRAN lines of their code and calculated the pressure on the obstacles in the medium through the momentum transferred by the reflected particles. Note, however, that one can use obstacles of any shape or size. One advantage of the CA methods is the ease with which such irregular systems with complex boundaries between the fluid and the solid can be studied. There have been a few papers in which CA methods have been used for simulating flow in various models of porous media. Balasubramanian et aI. (1987) used a triangular lattice in their CA simulations of flow in a porous medium in which a random fraction of sites of the lattice is blocked at random. The blocked sites represent the solid matrix of the medium. However, their study was limited to very high porosity (about 0.99). Rothman (1988) and Rothman and Keller (1988) also used CA models to simulate tlow in porous media employing a triangular lattice. In particular, Rothman (1988) used rectangular obstacles of various sizes to simulate the solid matrix of the medium. One interesting result of his work is that, in order that the CA limit, the simulation results approach the continuum mean size of the void area in the lattice must be at least twice the mean-free path of the particles (which is about 9 lattice bonds). This condition is obeyed in our simulations. This condition also implies that in order to simulate flow in low-porosity porous media, one has to use large lattices. Rothman (1988) did not study the dependence of the permeability on the porosity of the system. These authors also used relatively small lattices (of order of lo5 sites). The diffusivity calculations use the well-known result that the mean squared displacement R2 of a diffusing particle at time t, executing a random walk in the medium, is related to the diffusivity D by R2 = Dt apart from a lattice-dependent constant. In the physics literature, following de Gennes (1976), the diffusion problem (random walks in a disordered medium) is
of 30 circular
obstacles
in a channel.
2228
.
MUHAMMAD~AHIMI
and
often called the ant in a Inbyrinth, since the diffusers (the ants) try to find their way through the maze of various sites of the lattice. For our diffusion calculations we used a triangular or simple cubic lattice, both of which have coordination number 6. In the beginning of the simulations, a site is either selected with probability p to be open to diffusion, in which case it is called an allowed site and represents part of the pore space, or it is closed to diffusion, called a prohibited site, and represents a part of the solid matrix of the medium. This is perhaps the simplest method of incorporating a disordered morphology into the model. A diffusing ant takes a step along one of the six bonds that emanate from its present lattice site with a probability l/6 (or l/z for a network of coordination number z). If the newly-selected site is prohibited, then the ant stays at its present site; otherwise it moves to the selected neighbor site. However, more complex geometrical properties of the medium could, in principle, be incorporated into the simulations. For example, if a porous medium is characterized by its pore throat size distribution, then the throat sizes are assigned to the bonds between lattice sites. In this case the ant makes a transition to another allowed site with a probability proportional to the permeability of the bond between these two sites, where the permeability is proportional to the square of the pore size or effective channel radius (this law of transition probability can easily be derived by discretizing the continuum diffusion equation). Thus, anisotropy, stratification, or any other form of heterogeneity could be handled by this method and, in fact, in the present paper we simulate diffusion in a stratified porous medium. The diffusion calculations are also welldocumented in the literature (Pandey et al., 1984) for a Cyber 205 (or ETA 10) vector computer. The simplifications of the code of Pandey et al. for a Cray computer were made by Roman (1990) whose program we modified for our purposes. To explain the algorithm, we now give the innermost loop of the code. For simplicity, we consider the case when diffusion in Y-direction is monitored; the generalization to all directions is obvious. DO 50 II = 1, IANTS IN(I1) = IFIX(6 * RANF( )) NEWPOS(I1) = IPOS(I1) + IW(IN(I1)) JW = NEWPOS(II)/64 IF((l _AND SHIFTR(MOCC(JW), NEWPOS(T1) - 64*JW)).EQ.O) GO TO 50 IPOS(I1) = NEWPOS(I1) IY(I1) = IY(I1) + NY(IN(I1)) 50 CONTINUE The first line gives the loop over the number IANTS of various diffusors or ants; typically we used about 1000 such ants. Since the various ants are diffusing independently of each other, this loop can be vectorized. The next line calls the built-in random number generator RANF and thus gives an index IN(I1) between 0 and 5, corresponding to one of the
D~ETRICH STALTFFER six lattice directions. NEWPOS is the ant’s new position; thus IW(IN( )) must give the shift in the coordinate due to a move in the direction given by the index IN. We numbered all sites consecutively from 0 to L3 - 1, and thus the array IW contains the elements 1, - 1, L, - L, L*, - L’, on the simple cubic lattice. For the triangular lattice the neighbor shifts in IW are 1, - 1, L, - L, L - 1 and - L + 1. Now we have to find the bit which tells us if the newly-selected site is allowed or prohibited for the diffusing ant. This information is stored in an array MOCC(JW), JW = 0, 1,2,. _ , L2/64 - 1. If, for example, the number NEWPOS of that new site is 150 and we store all sites in 64-bit computer words, the information is stored in the computer word MOCC(JW) with JW = 2, which describes sites 128 to 191. Thus site 150 is stored in bit number NEWPOS - 64* JW = 150 - 64*2 = 22, if the rightmost bit of a computer word has bit number 0 and the leftmost has bit number 63 (or 31 on computers with four-byte arithmetic). Shifting the appropriate computer word MOCC by 22 bits to the right (SHIFTRfunction) and reading off the then rightmost bit by a logical AND with the number 1, we find the status of the site in question. (While the computer may believe that MOCC is an integer, we use MOCC merely as a string of bits.) If that site is prohibited, we do nothing (GO TO 50); otherwise we shift the posltion IPOS of the ant to the new position NEWPOS. The distance IY in Y-direction of the ant from the origin of )) to its random walk is updated by adding NY(IN( it. where the six elements of the array NY have the values, 0, 0, 1, - l,O, and 0, depending on the six lattice directions in the simple cubic lattice. Similar tricks were used in the two-dimensional triangular lattice, where the distances IY are somewhat more complicated. Even if an ant starts in the center of the lattice, it may move towards its boundaries. We thus work with periodic (or helical) boundaries, i.e. we add or subtract L3 from the actual position TPOS if IPOS has become too small or too large, respectively. (We do not change the distance IY in such boundary corrections: thus IY can become much larger than L.) If we would check for such boundary effects after each time step we would loose speed; thus we added 21 buffer planes to both the top and the bottom of the lattice, containing the same information in MOCC as the 21 lowest and highest planes, respectively, of the lattice. Only at rare time intervals do we check for ants close to the boundary and change, if needed, their positions by L3. A 22nd buffer plane with only MOCC = 0 (prohibited sites) is added to both ends to deal with the rare events that between two boundary checks an ant has crossed the whole buffer area. Previous experience with various versions of this algorithm (see, Roman, 1990, for further literature) shows that effects from the boundaries can still be important even for lattices containing millions of sites. Thus storage in single bits is important to allow for large lattices at the expense of speed. Nevertheless,
Efficient simulation of flow and transport in porous media
a single Cray-YMP processor made nearly 4 million diffusion steps per second. The speed would be higher (Paetzold, 1990) had we stored only one site in each computer word. In the present paper, for diffusion simulations we used cubic lattices of sizes up to 10243, and triangular lattices of sizes mostly 9600 x 4500 sites. We also averaged the results over all ants and ten independent realizations of the lattices. Each diffusing ant was allowed to take up to 217 steps (where each step is equivalent to one unit of time). For flow simulations using the CA method, we used triangular lattices of sizes up to 2000 x 666. However, because of vectorization, the codes are very efficient, and larger lattices can be used, if necessary. A copy of the codes can be obtained from the authors.
3. RESULTS
AND
DISCUSSIONS
We simulated diffusion in a stratified porous medium, using two- and three-dimensional lattices as models of the pore space. The motivation for this study was partly the works of Matheron and de Marsily (1980) and Bouchaud et al. (1990) who investigated random watk of a tracer particle in a stratified medium in which there are spatially inhomogeneous, but correlated, velocity fields. This is an important problem and is relevant to many fluid flow phenomena in porous media, such as oil recovery processes and ground water transport in geological aquifers (see, e.g., Dagan, 1987, and references therein). In such systems, the steady-state velocities in the longitudinal X-direction (parallel to the strata) generate a net bias in the motion of the particle. Although the average of this bias over an infinite number of strata is zero, it can take on nonzero and fluctuating values if the number of the strata is finite. Matheron and de Marsily (1980) and Bouchaud et al. (1990) showed that, if the number of strata is infinite, then, the mean-squared displacement of the random walker increases with time as R2 - t312. That is, one has superd5fSusion since the growth of R* with I is faster than linear. This superdiffusion behavior seems to be a generic property of random walks in a random velocity field in porous media with strong spatial correlations (Sahimi, 1987; Sahimi and Imdakm, 1988). As a result of this superdiffusive behavior, the diffusion coefficient and, thus, the dispersivity which is the ratio of the diffusivity and macroscopic velocity, increase with time (and distance), in agreement with field observations (Dagan, 1987). We thus studied the analogous problem of diffusion in a stratified medium with no velocity field, that is, the diffusers move in a static medium. We used the triangular and simple-cubic lattices, and divided them into N parallel strata of equal thickness W. It would not be difficult to assign a different thickness to each stratum, if experimental information is available. The planes of the strata were CES 45:9-E
2229
parallel to the X-direction. We studied a simple heterogeneous system in which, to assign random properties to the strata, we designated at random a fraction p of the sites of each stratum as the allowed or open sites, and blocked the rest. For the different strata, we distributed p uniformly between pc, the percolation threshold of the lattice (Stauffer, 1985), and unity, where pc = l/2 and 0.3116, for the triangular and cubic lattices, respectively. In this way, each stratum possesses its own effective properties and, as discussed below, over much of the range of p, the effective diffusivity of each stratum is proportional to p. The model porous medium generated in this way is macroscopically heterogeneous and anisotropic. Thus, we measured the effective diffusivities in both the X- and Y-directions (perpendicular to the strata). We also varied the number (and thus the thickness) of the strata in order to study its effect on the effective diffusivities. We emphasize again that, if appropriate experimental data are available, it would not be difficult at all to incorporate other kinds of heterogeneities into the model, and assign the effective properties The of the strata according to a given distribution. most time consuming part of the code is the execution of the random walk, and the generation of the lattice with a given spatial distribution of heterogeneities consumes only a small fraction of the total execution time. Figure 3 represents the variations with time of the effective diffusivities D, and D, in the triangular lattice, while Fig. 4 displays the same quantities in the cubic network. As can be seen, after a long enough time, D, reaches an essentially constant value, and no superdiffusion is observed. In contrast with D,, D, appears to decrease with time, never reaching a constant value. This behavior can be explained as follows. In the Y-direction, the strata effectively act as conductors in series. Since p is uniformly distributed between pc and 1, for a large enough system, at least one of the strata will have a value p = pf, which means that its effective diffusivity (or conductivity) will essentially be zero (since diffusion coefficients vanish at the
DB; 0.6 -
‘:
r:r
8::.-_-
n
x x
0.4-
x
Y
x xlY
0.2 -
I 0
1
I 2
I 3
I 4
I 5
Time
Fig. 3. Plot of diffusivity (normalized to unity for p = 1) vs decadic logarithm of times for the triangular network. The round symbols refer to D,, the crosses to D,.
2230
SAHIMI and DIETRICH STAUFFER
MUHAMMAD
. * . x
. a
02
0
c
. I
x x
I 1
xx
I
I
I
I
2
3
4
5
6
Time
Fig. 4. Plot of diffusivity (normalized to unity for p = I) vs decadic logarithm or time, for Ihc cubic network. The round symbols refer to D,, the crosses to D,.
percolation threshold of an infinite system). Therefore, after a long enough time, the overall D, should be zero. Had we distributed p in the range from a to 1, where a is some constant larger than pC, or not assigned the effective properties of each stratum according to a percolation algorithm, D, would also have reached a constant value after a long enough time. Thus, the apparent variation of D, with time is nothing but an indication that D, will eventually vanish. This can also be shown by a more formal and mathematical argument, but we believe the above argument is clear enough. However, there are certain situations in which D, and D, may never attain a constant value, even if the effective properties of the strata are not assigned according to a percolation algorithm used here. For example, if the effective permeabilities of the strata are correlated with a very large correlation length r, then, the effective diffusivities will never attain constant values, because the root meansquared displacement R of the ants wilI at most be of the order of the permeabiIity correlation length. Diffusive behavior will be reached only if R $ c. The advantage of the method we have described here is that such complications can be studied efficiently. In Fig. 5 we present the variations of the bdk effective diffusivities (i.e., when the entire system is made ofjust one stratum) with p in both the triangular and cubic lattices. The diffusivities have been normalized with respect to their values at p = 1. In this case of course, D, = D, = D. It is clear that, except for a small region near pe, D varies essentially linearly with p, and can be well represented by D = 2p ~ 1 in the triangular lattice, and by D = (3~ - I)/2 in the cubic lattice, which are similar to the efTectivemedium approximation (see, e.g., Kirkpatrick, 1973; Joy and Strieder, 1978). These formulae can be used to estimate D, for both the triangular and cubic lattices. Since in the X-direction the strata effectively act as conductors in parallel, the overall D, is roughly the average of D,(p) for each stratum, where the averaging is taken over all values of p (distributed uni-
Fig. 5. Bulk diffusivity D in two (x ) and three (dots) dimensions. The results are the averages of ten 1024 x 1024 x 1024 lattices
08-
0.6-
a
.
.
.
.
.
. 0.4-
02-
0
I IO
1 20
I 30
I 40
I 50
I 60
I 70
I3
N
Fig. 6. The diffusivity D, as a function of the number Iv of strata; averaged over ten 384 x 384 x 150 cubic lattices.
formly between zero and one). Using the above approximations for D,, we then estimate that D, = 0.5 for both two- and three-dimensional lattices. These estimates are valid if the strata are infinitely thick. However, in strata with a finite thickness the effect of the interfaces between the strata can be strong, so that D, cannot be estimated precisely by the sort of argument that we give here. Figure 6 shows the variations of D, with the number of strata in the cubic lattice. When the number of strata is very small, so that the effect of the interfaces between various strata is small, the porous medium behaves similar to its bulk behavior discussed above. Figure 6 shows that D, is roughly 0.5, in agreement with our estimate given above. On the other hand, if we consider the opposite limit in which the thickness of each stratum is exactly one (i.e., the total number of strata is the same as the thickness of the medium), the diffusing ants “see” the system approximately as if it has no strataand is totally random, since the time that they spend for travelling between the strata is large compared to the time they spend in each stratum.
Efficient simuIation of flow and transport in po~oub media Thus, in this limit the effective value of D, should also be close to the average value of about 0.5. Figure 5 shows that D, Y 0.52, in reasonable agreement with our estimate. Of course, we do not expect DX to be exactly the same as the average value estimated above with the aid of the effective medium approximation, if the strata are not too thin or too thick, since the interfaces between various strata do influence D,. In between these two limits, D, takes on values that are roughly constant, as shown in Fig. 6. We found roughly the same behavior for the triangular lattice and, thus, the results are not shown. We also simulated flow in a two-dimensional model porous medium using the CA method described above. We considered a system in which circular obstacles of the same radius were arranged in an essentially periodic way in a triangular lattice, and flow takes place through the channels between the circles (see Fig. 2). As such the phenomenon we studied is the two-dimensional analog of the problem of flow between periodic arrays of spheres that has been studied by many authors (see, e.g., Zick and Homsy, 1982, and references therein). We used three rows of circles in each of which ten circles of equal radius Y were placed. We varied r (and thus the porosity) in order to study its effect on the permeability k of the system. The results were obtained with a 2000 x 666 lattice. As r increases, the permeability k decreases and eventually vanishes as the circles start to overlap (i.e., as the system reaches its percolation threshold). We also used triangular lattices of size 900 x 300, and the results were in qualitative agreement with those of the larger system, though with lower permeabilities. In Fig. 7 we compare the permeability of our system with that calculated by Brosa and Stauffer (1991). These authors used the same lattice of the same size as ours, but distributed the circles in the system at random so that they can, in principle, overlap. As shown, at high values of porosities there is very little differ-
8
9
7
.
6 Y
5
2231
ence between the two permeabilities. However, as porosity decreases, the difference increases. For low values of porosity, the two permeabilities do not resemble one another at all. This figure indicates that the spatially periodic systems that have been used by many to model disordered porous media (for a review see, Nitsche and Brenner, 1989) is not an appropriate model of disordered porous media, especially for intermediate and low values of porosity. Indeed, in the case of a spatially periodic system, one has a predictable configuration of the pore space, whereas in the case of a random porous medium one may have large pores or channels through which a lot of fluid can flow, or large connected clusters of overlapping obstacles that allow little fluid to flow, giving rise to a percolation phenomenon. The spatially periodic models cannot predict the existence of a nontrivial percolation threshold (pc # 1 or 0) for the effective porosity of the medium. Before closing this section, we should point out that a long standing and unsolved problem is the exact relation between the electrical conductivity of a Auidsaturated porous medium and permeability of the same system. Of course, according to Einstein’s relation, the electrical conductivity is proportional to the diffusivity of the system (that is calculated here). Although various empirical relations between the two properties have been suggested, none has been successful over the entire range of porosities (see, e.g., Johnson et al., 1987, and references therein). In the context of the problems treated here, it suffices to say that, while the diffusivity of the system takes on a finite value at unit porosity, the permeability of the system is very large at this point. Although both quantities vanish at the percolation threshold of the system, there appears to be no relation between the two quantities. Johnson er al. (1987) have proposed that k = 11*/8F, where F is the formation factor of the medium, which is directly proportional to the diffusivity, and L is a well-defined parameter which is related to the dynamically-connected pores of the medium, and is roughly the ratio of the total volume and surface area of such pores. The CA method presented here can be used to calculate 1 for any configuration of the pore space and, therefore, investigate the relation between k and D. This matter will be studied in a future paper.
.+
.
4
.
:I+ , ,... :;....fy=; oq65
0.7
0.75
0.8
0.85
0.9
0.95
1
P
Fig. 7. Permeability k (crosses, arbitrary units) vs porosity p for 30 circular obstacles in three periodic layers, as shown in Fig. 2. The dots give the empirical law, k N (p - 2/3)/ (1 - p). found by Brosa and Stauffer (1990) for p above 3/4 for the same system but with the circles distributed at random locations.
4. SUMMARY
AND
CONCLUSIONS
In this paper, we have described two efficient computational algorithms for studying flow and transport processes in porous media, using lattice models of pore space. The two methods are ideally suited for use with a vector computer, and have made it possible to study flow and transport processes in porous media in great details. Although the two methods are not directly related to one another, the random-walk method can be used as a test of the accuracy of CA methods for simulating diffusion processes in porous media, since CA can be used for both ditTusion and
2232
MUHAMMADSAHIMI
convection.
This
comparison
and DIETRICH STAUFFER
will be the subject of
a future paper. The random-walk
t
method that we described here for studying transport processes in more complex systems. For example, to simulate the effect of a velocity field and study the problem that was investigated by Matheron and de Marsily (1980) and Bouchaud et al. (1990), one may introduce a bias in the random walk of the ants such that the probability of movement in the positive Xdirection is larger than those for the other directions. One may also use the algorithm to study transport in a porous medium with double porosity, i.e. a medium with large fractures and much smaller pores. To study such problems, we only need to introduce the proper configuration of the pore space for the ants to diffuse in and this, as discussed above, consumes only a small fraction of the total execution time. Thus, it is possible to study very large systems and investigate such phenomena in as much detail as needed. This will be the subject of a future paper. Elsewhere we have used the CA method to calculate the dynamic permeability of a porous medium (Duarte and Sahimi, 1991).
can be appropriately
Finally, we should emphasize again that one of the
Acknowledgefnents-The work of D.S. was partially supported by BMFT grant 0326657D. The work of M.S. was partly supported by the USC Center for Oil Recovery Research, Monterey Formation project. The contribution of the donor companies Elf Aquitaine, Oryx, Phillips Petroleum, Shell Development, Unocal and Texaco is gratefully acknowledged. We are also grateful to the San Diego Supercomputer Center for computer time. One of us (MS.) would like to thank Supercomputer Center HLRZ-KFA for warm hospitality during Summer of 1990. NOTATION DX DY F
k N P P. r
W
displacement time thickness of a stratum
modified
main advantages of CA methods for simulating flow in porous media is the ease with which they can be applied to complex geometries. Another chief virtue of CA methods is their adherence to the Navier-Stokes equations, thus allowing the full range of Bow behavior in porous media to be studied. In a future paper, we shall compare the predictions of CA methods for permeability of model porous media with various lower and upper bounds that have been derived for k (see Torquato, 1987 for a review) in order to assess the accuracy of such bounds. Since the completion of this work, Kohring (1991) has improved the efficiency of our computer program for CA method by a factor of about ten. Thus, simulation of flow and transport in porous media can be done even murt: eflkienlly than that discussed here.
D
R
bulk diffusivity diffusivity in the X-direction diffusivity in the Y-direction formation factor permeability number of strata porosity percolation threshold obstacle radius
Greek letters R parameter
5
pores correlation
related
to
dynamically
length of permeability
connected
heterogeneities
REFERENCES Balasubramanian, IL, Hayot, F. and Saam, W. F., 1987, Darcy’s law from lattice-gas hydrodynamics. Phys. Rev. A 36, 2248-2253. Barber, M. N. and Ninham, B. W., 1970, Random and Restricted Walks: Theory and Applications. Gordon and Breach, New York. Binder, K. (ed.), 1987, Applications of the Monte Carlo Method in Sratistical Physics. Springer, Berlin. Bouchaud, J.-P., Georges, A., Koplik, J., Provata, A. and Redner, S., 1990, Superdiffusion in random velocity fields. Phys. Rev. L&t. 64, 2503-2506. Brosa, U. and Stauffer, D., 1989, Vectorized multisite coding for hydrodynamic cellular automata. J. Stat. Phys. 57, 399-403. Brosa, U. and Stauffer, D., 1991, Simulation of flow through a 2D random porous medium. J. Stat. Phys. 62. Dagan, G., 1987, Theory of solute transport by groundwater. Annu. Rev. Fluid Me&. 19, 183-215. d’Humieres, D. and Lallemand, P., 1986, Lattice gas automata for fluid mechanics. Physica 14OA, 326-335. d’Humieres, D., Lallemand, P. and Frisch, 1986, Lattice gas models for 3D hydrodynamics. Europhys. Left. 2,291-297. Doolen, G. (ed.), 1991, Proceedings of the NATO Advanced Research Workshop on Lattice Gas Methods for PDE’s: Theory, Application and Hardware. Physica D 47. Duarte, J. A. M. S. and Brosa, U., 1990, Viscous drag by cellular automata. J. Stat. Phys. 59, 501-508. Duarte, J. A. M. S. and Sahimi, M., 1991, Dynamic permeability of a porous medium by celluIar automata. Preprint. Frisch, U., Hasslacher, B. and Pomeau, Y., 1986, Lattice-gas automata for Navier-Stokes eauations. Phvs. Reo. Lett. 56, 1505-1.508. de Gennes. P. G., 1976, La percolation: un concept unificateur. La Recherche 7, 914-927. Jerauld, G. R., Striven, L. E. and Davis, H. T., 1984, Percolation and conduction on the 3D Voronoi and regular networks: a second case study in topological disorder. J. Phys. C 17, 3429-3439. Johnson, D. L., Koplik, J. and Dashen, R., 1987, Theory of dynamic permeability and tortuosity in fluid-saturated porous media. J. Fluid Mech. 176, 379-402. Joy, T. and Strieder, W:, 1978, Effective medium theory of percolation in a random site triangular conductance network. J. Phys. C 11, L867-L870. Kirkpatrick, S., 1973, Percolation and conduction. Rev. Mod. Phys. 45, 574-588. Kuhring, G. A., 1991, Calculation of the pcmeability of porous media using hydrodynamic cellular automata. J. Stat. Phys. 62. Mason, E. A. and Malinauskas, A. P., 1983, Gas Transport in Porous Media: the Dusty-gas Model. Elsevier, Amsterdam. Matheron, G. and de Marsily, G., 1980, Is transport in porous media always diffusive? A counterexample. Water Resow. Res. 16, 901-917. Nitsche, L. C. and Brenner, H., 1989, Eulerian kinematics of flow through spatially periodic models of porous media. Archs. Ration. Mech. Analysis 107, 2255292. Paetzold, O., 1990, Diffusion of lattice gases without double
Efficient
simulation
of flow and transport
occupancy on 3-D percolation lattices. J. Stat. Phys. 61, 495-500. Pandey, R. B., Stauffer, D., Margolina, A. and Zabolitzky, J. G., 1984, Difftision on random systems above, below and at their percolation threshold in two and three dimensions. J. Stat. Phys. 34, 427-450. Roman, H. E., 1990, Diffusion in three-dimensional random systems at their percolation threshold. J. Stat. Phys. 58, 375-382. Rothman, D. H., 1988, Cellular-automata fluids: a model for flow in porous media. Geophysics 53,509-518. Rothman, D. H. and Keller, J. M., 1988, Immiscible cellularautomata fluids. J. Stat. Phys. 52, 1119-1127. Sahimi, M., 1987, Hydrodynamic dispersion near the percolation threshold: scaling and probability densities. .7. Phys. A 20, L1293-L1298. Sahimi, M., Gavalas, G. R. and Tsotsis, T. T., 1990, Statistical and continuum models of fluid-solid reactions in porous media. Chem. Engng Sci. 45, 1443-1502.
in porous
media
2233
Sahimi, M. and Imdakm, A. O., 1988, The effect of morphological disorder on hydrodynamic dispersion in flow through porous media. J. Phys. A 21, 3833-3870. Scheidegger, A. E., 1974, The Physics ofFluid Through Porous Media, 3rd Edition. University of Toronto Press, Toronto. Stauffer, D., 1985, Introduction to Percolation Theory. Taylor and Francis, London. Torquato, A., 1987, Thermal conductivity of disordered heterogeneous media from the microstructure. Rev. Chem. Engng 4, 15 l-204. Whitaker, S., 1986, Flow in porous media: the governing equations for immiscible, two phase flow. Transport in Porous Media 1, 105-125. Wolfram, S., 1986, Theory and AppIications of Cellular Automata. World Scientific, Singapore. Zick, A. A. and Homsy, G. M., 1982, Stokes flow through periodic arrays of spheres. J. Fluid Mech. 115, 13-26.
EFFICIENT
am!-2509/91 S3.03+ 0.00 0 1991 Pergamon Press plc
SIMULATION OF POROUS
FLOW AND MEDIA
TRANSPORT
IN
MUHAMMAD SAHIMI + and DIETRICH STAUFFER* Supercomputer Center HLRZ, c/o KFA, D-5170 Jiilich 1, Germany (First receiued
in revised form 20 November 1990)
10 Auyusf 1990; accepted
Abstract-We describe two efficient computer simulation methods for investigating flow and transport processes in porous media. The first method, which used random-walk techniques, is used for estimating effective diffusion coefficients in disordered porous media. The second method, which is based on lattice gases and celIular automata models, is used for calculating the effective permeability of a porous medium. We describe in some detail the computer algorithms that employ these methods. Etticlent use of vectonzation enables us to simulate diffusion in disordered porous media, representedby two- or three-dimensional lattices, with up to IO9sites, and a speed of nearly four steps per microsecond and per Gray-YMP processor. Vectorization of the algorithm based on cellular automata also enables us to simulate flow in a twodimensional model porous medium with more than a million sites. As examples, we study dilksion and flow in stratified and macroscopically heterogeneous porous media and Row through periodic arrays of obstacles. The relation between the effectivediffusivity and permeability of such porous media is-also briefly
discussed.
1. INTRODUCTION The problems of flow and transport through media are of great theoretical and practical
porous
interest (see, e.g., Scheidegger, 1974; Mason and Malinauskas, 1983; and references therein). For example, in the petroleum industry, the reservoir engineers have been dealing with such problems since oil production began. Problems such as drainage and imbibition in soil, multiphase flow through trickle-bed reactors, mercury porosimetry for determining the pore size distribution of a porous catalyst, channeling in packed columns and ground water flow are but a few processes of interest to hydrologists, soil scientists, and chemical and petroleum engineers. Historically, flow and transport processes in porous media have been simulated starting from continuum equations in which the relevant variables are related to conserved quantities, such as local density, local momentum flux, etc. Macroscopic properties, such as the effective transport coefficients, are then defined as averages of the corresponding microscopic or local quantities (see, e.g., Whitaker, 1986; and references therein). The averages are taken over a volume which is small compared to the volume of the system, but large enough for the transport equations to hold when applied to the volume. At every point of a porous medium one uses the smallest such volume and, thereby, generates macroscopic field variables obeying equations such as Fick’s law of diffusion and Darcy’s law of flow in porous media. This approach has been useful for macroscopic modelling and has been used
‘Author to whom all correspondence should be addressed. Present and permanent address: Department of Chemical Engineering, University of Southern California, Los Angeles, CA 90089-1211, U.S.A. z Present address: Institute for Theoretical Physics, Cologne University, D-5CMl Cologne 41, Germany.
extensively in the past. However, this approach is purely phenomenological and provides little insight into how flow and transport processes depend on the morphology of the pore space, i.e. its topology or pore connectivity, and geometry which is usually represented by a pore size distribution. In fact, in order to use such averaged transport equations, transport coefficients must a priori be known. Therefore, if the morphology of the pore space is strongly chaotic or disordered, a continuum approach may not be a quantitative method of studying transport processes in porous media. But when applicable, one usually uses numerical simulations to solve the continuum equations and obtains a few of the quantities of interest. Such numerical methods often involve some kind of discretization of the continuum equations. The second method of studying flow and transport processes in porous media is based on network or lattice models of porous media. One models the medium as a network of bonds, which represent the pore throats of the actual medium, and sites, which represent the pore bodies of the real system. The chaotic morphology of the medium can be incorporated into the network model by assigning randomly-selected effective sizes to the bonds and sites, to mimic the disordered geometry of the medium, and by using networks of random local coordination number, i.e. the number of bonds connected to the same site, to represent the chaotic topology of the pore space. It is often true that, if the coordination number of a regular network and the average coordination number of a random network are the same, the macroscopic properties of the two networks are very close (see, e.g., Jerauld et al., 1984). so that a network of fixed coordination number is often an adequate representation of the topology of the pore space. In a variation of the network models, one starts with a regular lattice and distributes (at random or in a prescribed fashion) 2225
2226
MUHAMMAD SAHIMI and DIETRICH STAUFFER
inclusions in the system in the form of disks, spheres, ellipsoids, or any other desired shape. The effective size of the inclusions can be fixed or a distributed quantity, and flow and transport take place through the channels between the overlapping or nonoverlapping inclusions. This is the model that we use here to study flow in porous media. Sahimi et al. (1990) have recently reviewed such network modelling approaches to the problems of estimating the effective transport properties of porous media. It can also be shown that there is a close analogy between such network or lattice models and the set of discretized equations that one obtains in the numerical simulations of the continuum equations, so that the network models can be justified on this ground also. Recent advances in computer simulation techniques (see, e.g., Binder, 1987, especially Chapter 8) have now made the network or lattice models an efficient tool of investigating flow and transport processes in porous media. In this paper, we describe two highly efficient methods for studying such processes in porous media. The first method is based on randomwalk techniques (see, e.g., Barber and Ninham, 1970) by which one can obtain reliable estimates of diffusivities and conductivities of disordered porous media. The second method is based on cellular automata (CA) which, in the present context, are essentially the discrete solution of the Navier-Stokes equations. We use the CA method for estimating the effective (Darcy) permeability of a porous medium. The random-walk methods of estimating the transport properties are of course well known. However, recent advances in supercomputer technology have made it possible to turn these methods into efficient computational tools, and study flow and transport processes in lattice or network models of porous media which can have a realistic (and often complex) geometry, using very large systems with millions of sites and bonds. The plan of this paper is as follows. In the next section we describe in some detail the random-walk and CA approaches, and we discuss how these methods can be vectorized for use in a vector computer. In section three, we use these methods to study flow and transport processes in stratified porous media, a highly important problem (see, e.g., Matheron and de Marsily, 1980, Bouchaud et al., 1990), and flow through periodic arrays of obstacles, a long-standing problem of considerable interest (see, e.g., Zick and Homsy, 1982). The paper is summarized in Section 4, where we also discuss the application of these methods to more complex problems. 2.
SIMULATION
simulations on vector or parallel computers like a Cray. A primitive example of CA is Prussian “bureaucracy”, where each site at the next time step assumes the state its superior site has at the present time step. For hydrodynamics, such automata are too simple, and actual computer codes for twodimensional systems use 128 or 256 different states for each site, i.e. 7 or 8 computer bits. The CA approximation of the Navier-Stokes equation in two dimensions (Frisch et al., 1986; Doolen, 1991) is based on particles either resting on the sites of a large triangular lattice, or moving with unit velocity on one of the six bonds emanating from each lattice site. If particles hit each other they are scattered according to the laws of momentum conservation. For example, one particle hitting another particle which is at rest, may be scattered such that its direction changes by MY’, and that the previously resting particle moves in a direction inclined at - 60” with respect to the direction of the incoming velocity (Fig. 1). The inverse process is also possible. If a particle hits a solid wall, it is reflected by 180” (no slip boundary condition). The triangular lattice geometry ensures that in this momentum transfer the total momentum vector is conserved. It also approximates the Navier-Stokes equations better than the square lattice in the continuum limit, because of the isotropy of the fourth-order tensor arising in this limiting process. The extension to three dimensions is more complex since no regular three-dimensional lattice is isotropic and, thus, in the continuum limit, one has spurious terms, in addition to those in the Navier-Stokes equation, which are caused by the anisotropy of the lattice. However, there are now several methods of circumventing this difficulty. For example, one can either use disordered neta three-dimensional and topologically work, such as the Voronoi network which is macroscopically isotropic. Alternatively, one can use (d’Humieres et al., 1986) a four-dimensional facecentered hypercubic lattice (FCHC). For this lattice, which has a coordination number of 24, all pairwise symmetric fourth-order tensors are isotropic and, thus, one is in a position to simulate four-dimensional Navier-Stokes equation on such a lattice. One observes that any solution of the four-dimensional equation which does not depend on the fourth dimension is a solution of the three-dimensional equations. This suggests the use of a FCHC lattice which wraps around periodically in the fourth direction. One can actually use a lattice which is only one lattice unit
high-speed
METHODS
Cellular automata (see, Wolfram, 1986, for a collection of reviews) are large lattices where each lattice site can be in one of several discrete states. The state of each site at the next time step is determined completely by the present state of the neighboring sites. Thus both time and space are discrete, and connections are between neighbors only, ideal conditions for
Fig. 1. Example for a transition rule in the cellular automata method. A particle moving with unit velocity hits a particle at rest, and both are scattered such that the total momentum is conserved in the triangular lattice.
Efficient
simulation of Row and transport
long in the fourth dimension and, therefore, have an effectively three-dimensional spatial structure. The disadvantage of this method is that, although the fourth dimension is very thin, the discrete velocities still have components in all directions; therefore, the model is somewhat bit intensive (24 or 25 bits per site). A third approach (d’Humieres et al., 1986) is to use a three-dimensional cubic lattice in which the particles move with speed zero, one and $ (instead of unit velocity in the two-dimensional lattices). This model uses only 19 bits per site. Despite the discrete nature of CA in both two and three dimensions, these models are capable of exhibiting rich macroscopic complexity such as turbulence (d’Humieres and Lallemand, 1986). In the present paper, the porous mednnn is put into a channel with fixed walls. Obstacles are then distributed in between the walls. In the limiting case of unit porosity (with no obstacles) a simple Poiseuille Aow is expected at low enough Reynolds numbers. At the beginning of the simulation, we construct a transition table which tells us how one of the present states (determined by velocity of the incoming particles) is transformed into the next state, determined by the outgoing velocities. This table contains all of the possible states (which, in our simulations, number 2’ = 256). We then start with a biased random selection of molecular directions, where the bias is such that the average local velocities obey the Poiseuille profile, and simulate the scattering of the molecules by the obstacles of which the porous medium consists of, using the transition table we constructed. There can be two-, three- and four-body collisions at a site, which conserve energy and momentum. The rules are such that no new particle is created, and in one time step all particles on the lattice move first to different site, and then undergo collisions at the new site. One iteration consists of updating all lattice sites according to the transition table, and of the order of lo3 iterations are usually needed to reduce the average velocity to a value which no longer diminishes significantly in further iterations. In this paper, data were taken using 3000 such iterations. The CA method can be used for any porosity, but for low values of porosity, one needs a larger number of iterations to achieve stable results. In these simulations, we arrange circular obstacles in regular periodic arrays as shown in Fig. 2. Fluid flow then takes place in between nonoverlapping
Fig. 2. Flow between
a periodic
arrangement
in porous
media
2227
circles, and from the pressure and velocity distributions the effective permeability of the medium is calculated, just as in the case of one single circle or many randomly-placed and partially overlapping circles (Brosa and Stauffer, 1990; Duarte and Brosa, 1990). Various computational implementations of these ideas are discussed in the literature (Doolen, 1991); we follow Brosa and Staufier (1989) who published the innermost FORTRAN lines of their code and calculated the pressure on the obstacles in the medium through the momentum transferred by the reflected particles. Note, however, that one can use obstacles of any shape or size. One advantage of the CA methods is the ease with which such irregular systems with complex boundaries between the fluid and the solid can be studied. There have been a few papers in which CA methods have been used for simulating flow in various models of porous media. Balasubramanian et aI. (1987) used a triangular lattice in their CA simulations of flow in a porous medium in which a random fraction of sites of the lattice is blocked at random. The blocked sites represent the solid matrix of the medium. However, their study was limited to very high porosity (about 0.99). Rothman (1988) and Rothman and Keller (1988) also used CA models to simulate tlow in porous media employing a triangular lattice. In particular, Rothman (1988) used rectangular obstacles of various sizes to simulate the solid matrix of the medium. One interesting result of his work is that, in order that the CA limit, the simulation results approach the continuum mean size of the void area in the lattice must be at least twice the mean-free path of the particles (which is about 9 lattice bonds). This condition is obeyed in our simulations. This condition also implies that in order to simulate flow in low-porosity porous media, one has to use large lattices. Rothman (1988) did not study the dependence of the permeability on the porosity of the system. These authors also used relatively small lattices (of order of lo5 sites). The diffusivity calculations use the well-known result that the mean squared displacement R2 of a diffusing particle at time t, executing a random walk in the medium, is related to the diffusivity D by R2 = Dt apart from a lattice-dependent constant. In the physics literature, following de Gennes (1976), the diffusion problem (random walks in a disordered medium) is
of 30 circular
obstacles
in a channel.
2228
.
MUHAMMAD~AHIMI
and
often called the ant in a Inbyrinth, since the diffusers (the ants) try to find their way through the maze of various sites of the lattice. For our diffusion calculations we used a triangular or simple cubic lattice, both of which have coordination number 6. In the beginning of the simulations, a site is either selected with probability p to be open to diffusion, in which case it is called an allowed site and represents part of the pore space, or it is closed to diffusion, called a prohibited site, and represents a part of the solid matrix of the medium. This is perhaps the simplest method of incorporating a disordered morphology into the model. A diffusing ant takes a step along one of the six bonds that emanate from its present lattice site with a probability l/6 (or l/z for a network of coordination number z). If the newly-selected site is prohibited, then the ant stays at its present site; otherwise it moves to the selected neighbor site. However, more complex geometrical properties of the medium could, in principle, be incorporated into the simulations. For example, if a porous medium is characterized by its pore throat size distribution, then the throat sizes are assigned to the bonds between lattice sites. In this case the ant makes a transition to another allowed site with a probability proportional to the permeability of the bond between these two sites, where the permeability is proportional to the square of the pore size or effective channel radius (this law of transition probability can easily be derived by discretizing the continuum diffusion equation). Thus, anisotropy, stratification, or any other form of heterogeneity could be handled by this method and, in fact, in the present paper we simulate diffusion in a stratified porous medium. The diffusion calculations are also welldocumented in the literature (Pandey et al., 1984) for a Cyber 205 (or ETA 10) vector computer. The simplifications of the code of Pandey et al. for a Cray computer were made by Roman (1990) whose program we modified for our purposes. To explain the algorithm, we now give the innermost loop of the code. For simplicity, we consider the case when diffusion in Y-direction is monitored; the generalization to all directions is obvious. DO 50 II = 1, IANTS IN(I1) = IFIX(6 * RANF( )) NEWPOS(I1) = IPOS(I1) + IW(IN(I1)) JW = NEWPOS(II)/64 IF((l _AND SHIFTR(MOCC(JW), NEWPOS(T1) - 64*JW)).EQ.O) GO TO 50 IPOS(I1) = NEWPOS(I1) IY(I1) = IY(I1) + NY(IN(I1)) 50 CONTINUE The first line gives the loop over the number IANTS of various diffusors or ants; typically we used about 1000 such ants. Since the various ants are diffusing independently of each other, this loop can be vectorized. The next line calls the built-in random number generator RANF and thus gives an index IN(I1) between 0 and 5, corresponding to one of the
D~ETRICH STALTFFER six lattice directions. NEWPOS is the ant’s new position; thus IW(IN( )) must give the shift in the coordinate due to a move in the direction given by the index IN. We numbered all sites consecutively from 0 to L3 - 1, and thus the array IW contains the elements 1, - 1, L, - L, L*, - L’, on the simple cubic lattice. For the triangular lattice the neighbor shifts in IW are 1, - 1, L, - L, L - 1 and - L + 1. Now we have to find the bit which tells us if the newly-selected site is allowed or prohibited for the diffusing ant. This information is stored in an array MOCC(JW), JW = 0, 1,2,. _ , L2/64 - 1. If, for example, the number NEWPOS of that new site is 150 and we store all sites in 64-bit computer words, the information is stored in the computer word MOCC(JW) with JW = 2, which describes sites 128 to 191. Thus site 150 is stored in bit number NEWPOS - 64* JW = 150 - 64*2 = 22, if the rightmost bit of a computer word has bit number 0 and the leftmost has bit number 63 (or 31 on computers with four-byte arithmetic). Shifting the appropriate computer word MOCC by 22 bits to the right (SHIFTRfunction) and reading off the then rightmost bit by a logical AND with the number 1, we find the status of the site in question. (While the computer may believe that MOCC is an integer, we use MOCC merely as a string of bits.) If that site is prohibited, we do nothing (GO TO 50); otherwise we shift the posltion IPOS of the ant to the new position NEWPOS. The distance IY in Y-direction of the ant from the origin of )) to its random walk is updated by adding NY(IN( it. where the six elements of the array NY have the values, 0, 0, 1, - l,O, and 0, depending on the six lattice directions in the simple cubic lattice. Similar tricks were used in the two-dimensional triangular lattice, where the distances IY are somewhat more complicated. Even if an ant starts in the center of the lattice, it may move towards its boundaries. We thus work with periodic (or helical) boundaries, i.e. we add or subtract L3 from the actual position TPOS if IPOS has become too small or too large, respectively. (We do not change the distance IY in such boundary corrections: thus IY can become much larger than L.) If we would check for such boundary effects after each time step we would loose speed; thus we added 21 buffer planes to both the top and the bottom of the lattice, containing the same information in MOCC as the 21 lowest and highest planes, respectively, of the lattice. Only at rare time intervals do we check for ants close to the boundary and change, if needed, their positions by L3. A 22nd buffer plane with only MOCC = 0 (prohibited sites) is added to both ends to deal with the rare events that between two boundary checks an ant has crossed the whole buffer area. Previous experience with various versions of this algorithm (see, Roman, 1990, for further literature) shows that effects from the boundaries can still be important even for lattices containing millions of sites. Thus storage in single bits is important to allow for large lattices at the expense of speed. Nevertheless,
Efficient simulation of flow and transport in porous media
a single Cray-YMP processor made nearly 4 million diffusion steps per second. The speed would be higher (Paetzold, 1990) had we stored only one site in each computer word. In the present paper, for diffusion simulations we used cubic lattices of sizes up to 10243, and triangular lattices of sizes mostly 9600 x 4500 sites. We also averaged the results over all ants and ten independent realizations of the lattices. Each diffusing ant was allowed to take up to 217 steps (where each step is equivalent to one unit of time). For flow simulations using the CA method, we used triangular lattices of sizes up to 2000 x 666. However, because of vectorization, the codes are very efficient, and larger lattices can be used, if necessary. A copy of the codes can be obtained from the authors.
3. RESULTS
AND
DISCUSSIONS
We simulated diffusion in a stratified porous medium, using two- and three-dimensional lattices as models of the pore space. The motivation for this study was partly the works of Matheron and de Marsily (1980) and Bouchaud et al. (1990) who investigated random watk of a tracer particle in a stratified medium in which there are spatially inhomogeneous, but correlated, velocity fields. This is an important problem and is relevant to many fluid flow phenomena in porous media, such as oil recovery processes and ground water transport in geological aquifers (see, e.g., Dagan, 1987, and references therein). In such systems, the steady-state velocities in the longitudinal X-direction (parallel to the strata) generate a net bias in the motion of the particle. Although the average of this bias over an infinite number of strata is zero, it can take on nonzero and fluctuating values if the number of the strata is finite. Matheron and de Marsily (1980) and Bouchaud et al. (1990) showed that, if the number of strata is infinite, then, the mean-squared displacement of the random walker increases with time as R2 - t312. That is, one has superd5fSusion since the growth of R* with I is faster than linear. This superdiffusion behavior seems to be a generic property of random walks in a random velocity field in porous media with strong spatial correlations (Sahimi, 1987; Sahimi and Imdakm, 1988). As a result of this superdiffusive behavior, the diffusion coefficient and, thus, the dispersivity which is the ratio of the diffusivity and macroscopic velocity, increase with time (and distance), in agreement with field observations (Dagan, 1987). We thus studied the analogous problem of diffusion in a stratified medium with no velocity field, that is, the diffusers move in a static medium. We used the triangular and simple-cubic lattices, and divided them into N parallel strata of equal thickness W. It would not be difficult to assign a different thickness to each stratum, if experimental information is available. The planes of the strata were CES 45:9-E
2229
parallel to the X-direction. We studied a simple heterogeneous system in which, to assign random properties to the strata, we designated at random a fraction p of the sites of each stratum as the allowed or open sites, and blocked the rest. For the different strata, we distributed p uniformly between pc, the percolation threshold of the lattice (Stauffer, 1985), and unity, where pc = l/2 and 0.3116, for the triangular and cubic lattices, respectively. In this way, each stratum possesses its own effective properties and, as discussed below, over much of the range of p, the effective diffusivity of each stratum is proportional to p. The model porous medium generated in this way is macroscopically heterogeneous and anisotropic. Thus, we measured the effective diffusivities in both the X- and Y-directions (perpendicular to the strata). We also varied the number (and thus the thickness) of the strata in order to study its effect on the effective diffusivities. We emphasize again that, if appropriate experimental data are available, it would not be difficult at all to incorporate other kinds of heterogeneities into the model, and assign the effective properties The of the strata according to a given distribution. most time consuming part of the code is the execution of the random walk, and the generation of the lattice with a given spatial distribution of heterogeneities consumes only a small fraction of the total execution time. Figure 3 represents the variations with time of the effective diffusivities D, and D, in the triangular lattice, while Fig. 4 displays the same quantities in the cubic network. As can be seen, after a long enough time, D, reaches an essentially constant value, and no superdiffusion is observed. In contrast with D,, D, appears to decrease with time, never reaching a constant value. This behavior can be explained as follows. In the Y-direction, the strata effectively act as conductors in series. Since p is uniformly distributed between pc and 1, for a large enough system, at least one of the strata will have a value p = pf, which means that its effective diffusivity (or conductivity) will essentially be zero (since diffusion coefficients vanish at the
DB; 0.6 -
‘:
r:r
8::.-_-
n
x x
0.4-
x
Y
x xlY
0.2 -
I 0
1
I 2
I 3
I 4
I 5
Time
Fig. 3. Plot of diffusivity (normalized to unity for p = 1) vs decadic logarithm of times for the triangular network. The round symbols refer to D,, the crosses to D,.
2230
SAHIMI and DIETRICH STAUFFER
MUHAMMAD
. * . x
. a
02
0
c
. I
x x
I 1
xx
I
I
I
I
2
3
4
5
6
Time
Fig. 4. Plot of diffusivity (normalized to unity for p = I) vs decadic logarithm or time, for Ihc cubic network. The round symbols refer to D,, the crosses to D,.
percolation threshold of an infinite system). Therefore, after a long enough time, the overall D, should be zero. Had we distributed p in the range from a to 1, where a is some constant larger than pC, or not assigned the effective properties of each stratum according to a percolation algorithm, D, would also have reached a constant value after a long enough time. Thus, the apparent variation of D, with time is nothing but an indication that D, will eventually vanish. This can also be shown by a more formal and mathematical argument, but we believe the above argument is clear enough. However, there are certain situations in which D, and D, may never attain a constant value, even if the effective properties of the strata are not assigned according to a percolation algorithm used here. For example, if the effective permeabilities of the strata are correlated with a very large correlation length r, then, the effective diffusivities will never attain constant values, because the root meansquared displacement R of the ants wilI at most be of the order of the permeabiIity correlation length. Diffusive behavior will be reached only if R $ c. The advantage of the method we have described here is that such complications can be studied efficiently. In Fig. 5 we present the variations of the bdk effective diffusivities (i.e., when the entire system is made ofjust one stratum) with p in both the triangular and cubic lattices. The diffusivities have been normalized with respect to their values at p = 1. In this case of course, D, = D, = D. It is clear that, except for a small region near pe, D varies essentially linearly with p, and can be well represented by D = 2p ~ 1 in the triangular lattice, and by D = (3~ - I)/2 in the cubic lattice, which are similar to the efTectivemedium approximation (see, e.g., Kirkpatrick, 1973; Joy and Strieder, 1978). These formulae can be used to estimate D, for both the triangular and cubic lattices. Since in the X-direction the strata effectively act as conductors in parallel, the overall D, is roughly the average of D,(p) for each stratum, where the averaging is taken over all values of p (distributed uni-
Fig. 5. Bulk diffusivity D in two (x ) and three (dots) dimensions. The results are the averages of ten 1024 x 1024 x 1024 lattices
08-
0.6-
a
.
.
.
.
.
. 0.4-
02-
0
I IO
1 20
I 30
I 40
I 50
I 60
I 70
I3
N
Fig. 6. The diffusivity D, as a function of the number Iv of strata; averaged over ten 384 x 384 x 150 cubic lattices.
formly between zero and one). Using the above approximations for D,, we then estimate that D, = 0.5 for both two- and three-dimensional lattices. These estimates are valid if the strata are infinitely thick. However, in strata with a finite thickness the effect of the interfaces between the strata can be strong, so that D, cannot be estimated precisely by the sort of argument that we give here. Figure 6 shows the variations of D, with the number of strata in the cubic lattice. When the number of strata is very small, so that the effect of the interfaces between various strata is small, the porous medium behaves similar to its bulk behavior discussed above. Figure 6 shows that D, is roughly 0.5, in agreement with our estimate given above. On the other hand, if we consider the opposite limit in which the thickness of each stratum is exactly one (i.e., the total number of strata is the same as the thickness of the medium), the diffusing ants “see” the system approximately as if it has no strataand is totally random, since the time that they spend for travelling between the strata is large compared to the time they spend in each stratum.
Efficient simuIation of flow and transport in po~oub media Thus, in this limit the effective value of D, should also be close to the average value of about 0.5. Figure 5 shows that D, Y 0.52, in reasonable agreement with our estimate. Of course, we do not expect DX to be exactly the same as the average value estimated above with the aid of the effective medium approximation, if the strata are not too thin or too thick, since the interfaces between various strata do influence D,. In between these two limits, D, takes on values that are roughly constant, as shown in Fig. 6. We found roughly the same behavior for the triangular lattice and, thus, the results are not shown. We also simulated flow in a two-dimensional model porous medium using the CA method described above. We considered a system in which circular obstacles of the same radius were arranged in an essentially periodic way in a triangular lattice, and flow takes place through the channels between the circles (see Fig. 2). As such the phenomenon we studied is the two-dimensional analog of the problem of flow between periodic arrays of spheres that has been studied by many authors (see, e.g., Zick and Homsy, 1982, and references therein). We used three rows of circles in each of which ten circles of equal radius Y were placed. We varied r (and thus the porosity) in order to study its effect on the permeability k of the system. The results were obtained with a 2000 x 666 lattice. As r increases, the permeability k decreases and eventually vanishes as the circles start to overlap (i.e., as the system reaches its percolation threshold). We also used triangular lattices of size 900 x 300, and the results were in qualitative agreement with those of the larger system, though with lower permeabilities. In Fig. 7 we compare the permeability of our system with that calculated by Brosa and Stauffer (1991). These authors used the same lattice of the same size as ours, but distributed the circles in the system at random so that they can, in principle, overlap. As shown, at high values of porosities there is very little differ-
8
9
7
.
6 Y
5
2231
ence between the two permeabilities. However, as porosity decreases, the difference increases. For low values of porosity, the two permeabilities do not resemble one another at all. This figure indicates that the spatially periodic systems that have been used by many to model disordered porous media (for a review see, Nitsche and Brenner, 1989) is not an appropriate model of disordered porous media, especially for intermediate and low values of porosity. Indeed, in the case of a spatially periodic system, one has a predictable configuration of the pore space, whereas in the case of a random porous medium one may have large pores or channels through which a lot of fluid can flow, or large connected clusters of overlapping obstacles that allow little fluid to flow, giving rise to a percolation phenomenon. The spatially periodic models cannot predict the existence of a nontrivial percolation threshold (pc # 1 or 0) for the effective porosity of the medium. Before closing this section, we should point out that a long standing and unsolved problem is the exact relation between the electrical conductivity of a Auidsaturated porous medium and permeability of the same system. Of course, according to Einstein’s relation, the electrical conductivity is proportional to the diffusivity of the system (that is calculated here). Although various empirical relations between the two properties have been suggested, none has been successful over the entire range of porosities (see, e.g., Johnson et al., 1987, and references therein). In the context of the problems treated here, it suffices to say that, while the diffusivity of the system takes on a finite value at unit porosity, the permeability of the system is very large at this point. Although both quantities vanish at the percolation threshold of the system, there appears to be no relation between the two quantities. Johnson er al. (1987) have proposed that k = 11*/8F, where F is the formation factor of the medium, which is directly proportional to the diffusivity, and L is a well-defined parameter which is related to the dynamically-connected pores of the medium, and is roughly the ratio of the total volume and surface area of such pores. The CA method presented here can be used to calculate 1 for any configuration of the pore space and, therefore, investigate the relation between k and D. This matter will be studied in a future paper.
.+
.
4
.
:I+ , ,... :;....fy=; oq65
0.7
0.75
0.8
0.85
0.9
0.95
1
P
Fig. 7. Permeability k (crosses, arbitrary units) vs porosity p for 30 circular obstacles in three periodic layers, as shown in Fig. 2. The dots give the empirical law, k N (p - 2/3)/ (1 - p). found by Brosa and Stauffer (1990) for p above 3/4 for the same system but with the circles distributed at random locations.
4. SUMMARY
AND
CONCLUSIONS
In this paper, we have described two efficient computational algorithms for studying flow and transport processes in porous media, using lattice models of pore space. The two methods are ideally suited for use with a vector computer, and have made it possible to study flow and transport processes in porous media in great details. Although the two methods are not directly related to one another, the random-walk method can be used as a test of the accuracy of CA methods for simulating diffusion processes in porous media, since CA can be used for both ditTusion and
2232
MUHAMMADSAHIMI
convection.
This
comparison
and DIETRICH STAUFFER
will be the subject of
a future paper. The random-walk
t
method that we described here for studying transport processes in more complex systems. For example, to simulate the effect of a velocity field and study the problem that was investigated by Matheron and de Marsily (1980) and Bouchaud et al. (1990), one may introduce a bias in the random walk of the ants such that the probability of movement in the positive Xdirection is larger than those for the other directions. One may also use the algorithm to study transport in a porous medium with double porosity, i.e. a medium with large fractures and much smaller pores. To study such problems, we only need to introduce the proper configuration of the pore space for the ants to diffuse in and this, as discussed above, consumes only a small fraction of the total execution time. Thus, it is possible to study very large systems and investigate such phenomena in as much detail as needed. This will be the subject of a future paper. Elsewhere we have used the CA method to calculate the dynamic permeability of a porous medium (Duarte and Sahimi, 1991).
can be appropriately
Finally, we should emphasize again that one of the
Acknowledgefnents-The work of D.S. was partially supported by BMFT grant 0326657D. The work of M.S. was partly supported by the USC Center for Oil Recovery Research, Monterey Formation project. The contribution of the donor companies Elf Aquitaine, Oryx, Phillips Petroleum, Shell Development, Unocal and Texaco is gratefully acknowledged. We are also grateful to the San Diego Supercomputer Center for computer time. One of us (MS.) would like to thank Supercomputer Center HLRZ-KFA for warm hospitality during Summer of 1990. NOTATION DX DY F
k N P P. r
W
displacement time thickness of a stratum
modified
main advantages of CA methods for simulating flow in porous media is the ease with which they can be applied to complex geometries. Another chief virtue of CA methods is their adherence to the Navier-Stokes equations, thus allowing the full range of Bow behavior in porous media to be studied. In a future paper, we shall compare the predictions of CA methods for permeability of model porous media with various lower and upper bounds that have been derived for k (see Torquato, 1987 for a review) in order to assess the accuracy of such bounds. Since the completion of this work, Kohring (1991) has improved the efficiency of our computer program for CA method by a factor of about ten. Thus, simulation of flow and transport in porous media can be done even murt: eflkienlly than that discussed here.
D
R
bulk diffusivity diffusivity in the X-direction diffusivity in the Y-direction formation factor permeability number of strata porosity percolation threshold obstacle radius
Greek letters R parameter
5
pores correlation
related
to
dynamically
length of permeability
connected
heterogeneities
REFERENCES Balasubramanian, IL, Hayot, F. and Saam, W. F., 1987, Darcy’s law from lattice-gas hydrodynamics. Phys. Rev. A 36, 2248-2253. Barber, M. N. and Ninham, B. W., 1970, Random and Restricted Walks: Theory and Applications. Gordon and Breach, New York. Binder, K. (ed.), 1987, Applications of the Monte Carlo Method in Sratistical Physics. Springer, Berlin. Bouchaud, J.-P., Georges, A., Koplik, J., Provata, A. and Redner, S., 1990, Superdiffusion in random velocity fields. Phys. Rev. L&t. 64, 2503-2506. Brosa, U. and Stauffer, D., 1989, Vectorized multisite coding for hydrodynamic cellular automata. J. Stat. Phys. 57, 399-403. Brosa, U. and Stauffer, D., 1991, Simulation of flow through a 2D random porous medium. J. Stat. Phys. 62. Dagan, G., 1987, Theory of solute transport by groundwater. Annu. Rev. Fluid Me&. 19, 183-215. d’Humieres, D. and Lallemand, P., 1986, Lattice gas automata for fluid mechanics. Physica 14OA, 326-335. d’Humieres, D., Lallemand, P. and Frisch, 1986, Lattice gas models for 3D hydrodynamics. Europhys. Left. 2,291-297. Doolen, G. (ed.), 1991, Proceedings of the NATO Advanced Research Workshop on Lattice Gas Methods for PDE’s: Theory, Application and Hardware. Physica D 47. Duarte, J. A. M. S. and Brosa, U., 1990, Viscous drag by cellular automata. J. Stat. Phys. 59, 501-508. Duarte, J. A. M. S. and Sahimi, M., 1991, Dynamic permeability of a porous medium by celluIar automata. Preprint. Frisch, U., Hasslacher, B. and Pomeau, Y., 1986, Lattice-gas automata for Navier-Stokes eauations. Phvs. Reo. Lett. 56, 1505-1.508. de Gennes. P. G., 1976, La percolation: un concept unificateur. La Recherche 7, 914-927. Jerauld, G. R., Striven, L. E. and Davis, H. T., 1984, Percolation and conduction on the 3D Voronoi and regular networks: a second case study in topological disorder. J. Phys. C 17, 3429-3439. Johnson, D. L., Koplik, J. and Dashen, R., 1987, Theory of dynamic permeability and tortuosity in fluid-saturated porous media. J. Fluid Mech. 176, 379-402. Joy, T. and Strieder, W:, 1978, Effective medium theory of percolation in a random site triangular conductance network. J. Phys. C 11, L867-L870. Kirkpatrick, S., 1973, Percolation and conduction. Rev. Mod. Phys. 45, 574-588. Kuhring, G. A., 1991, Calculation of the pcmeability of porous media using hydrodynamic cellular automata. J. Stat. Phys. 62. Mason, E. A. and Malinauskas, A. P., 1983, Gas Transport in Porous Media: the Dusty-gas Model. Elsevier, Amsterdam. Matheron, G. and de Marsily, G., 1980, Is transport in porous media always diffusive? A counterexample. Water Resow. Res. 16, 901-917. Nitsche, L. C. and Brenner, H., 1989, Eulerian kinematics of flow through spatially periodic models of porous media. Archs. Ration. Mech. Analysis 107, 2255292. Paetzold, O., 1990, Diffusion of lattice gases without double
Efficient
simulation
of flow and transport
occupancy on 3-D percolation lattices. J. Stat. Phys. 61, 495-500. Pandey, R. B., Stauffer, D., Margolina, A. and Zabolitzky, J. G., 1984, Difftision on random systems above, below and at their percolation threshold in two and three dimensions. J. Stat. Phys. 34, 427-450. Roman, H. E., 1990, Diffusion in three-dimensional random systems at their percolation threshold. J. Stat. Phys. 58, 375-382. Rothman, D. H., 1988, Cellular-automata fluids: a model for flow in porous media. Geophysics 53,509-518. Rothman, D. H. and Keller, J. M., 1988, Immiscible cellularautomata fluids. J. Stat. Phys. 52, 1119-1127. Sahimi, M., 1987, Hydrodynamic dispersion near the percolation threshold: scaling and probability densities. .7. Phys. A 20, L1293-L1298. Sahimi, M., Gavalas, G. R. and Tsotsis, T. T., 1990, Statistical and continuum models of fluid-solid reactions in porous media. Chem. Engng Sci. 45, 1443-1502.
in porous
media
2233
Sahimi, M. and Imdakm, A. O., 1988, The effect of morphological disorder on hydrodynamic dispersion in flow through porous media. J. Phys. A 21, 3833-3870. Scheidegger, A. E., 1974, The Physics ofFluid Through Porous Media, 3rd Edition. University of Toronto Press, Toronto. Stauffer, D., 1985, Introduction to Percolation Theory. Taylor and Francis, London. Torquato, A., 1987, Thermal conductivity of disordered heterogeneous media from the microstructure. Rev. Chem. Engng 4, 15 l-204. Whitaker, S., 1986, Flow in porous media: the governing equations for immiscible, two phase flow. Transport in Porous Media 1, 105-125. Wolfram, S., 1986, Theory and AppIications of Cellular Automata. World Scientific, Singapore. Zick, A. A. and Homsy, G. M., 1982, Stokes flow through periodic arrays of spheres. J. Fluid Mech. 115, 13-26.