Finite-size anisotropy in statistically uniform porous media
Anisotropy of the permeability tensor in statistically uniform porous media of sizes used in typical computer simulations is studied. Although such systems are assumed to be isotropic by default, we show that de facto their anisotropic permeability can give rise to significant changes of transport parameters such as permeability and tortuosity. The main parameter controlling the anisotropy is , being the ratio of the obstacle to system size. Distribution of the angle between the external force and the volumetric fluid stream is found to be approximately normal, and the standard deviation of is found to decay with the system size as , where is the space dimensionality. These properties can be used to estimate both anisotropy-related statistical errors in large-scale simulations and the size of the representative elementary volume.
A standard method of modeling a uniform, isotropic porous medium (e.g. a column of sand) is to place randomly many identical objects that are impermeable to fluid (e.g. solid spheres) in an initially empty volume Martys and Garboczi (1992); Koponen et al. (1997); Nabovati and Sousa (2007); Verberg and Ladd (1999); Kostek et al. (1992); Matyka et al. (2008); Başagaoglu et al. (2006); Zhang and Knackstedt (1995). Since the objects are placed uniformly in the whole system, one might expect that randomness in their exact locations is irrelevant in the sense that the bulk volumetric fluid stream will be parallel to the external force (e.g. gravitation). This would be the case if the system was large enough. However, in computer simulations and in artificial laboratory systems (used in particle image velocimetry measurements Zerai et al. (2005)), usually relatively small systems are utilized that contain at most a few thousands of ”grains”—far less than billions of sand grains in a typical experimental setup. Since randomly distributed grains tend to form channels of random orientations, small porous systems are very sensitive to local fluctuations of the grain distribution. Under such conditions the direction of the volumetric fluid stream can differ significantly from that of the external force. Consequently, a system that was supposed to be isotropic, may de facto be rather highly anisotropic. The aim of this paper is a detailed analysis of this phenomenon in a two-dimensional (2D) flow.
A porous medium is anisotropic to flow if the permeability tensor is anisotropic. Usually anisotropy of the permeability tensor is either assumed explicitly Bear (1972) or it appears naturally as an expected consequence of a microscopic model Scholes et al. (2007); Boudreau and Meysman (2006); Alam et al. (2006); Ohkubo (2008). In the former case one works entirely on a macroscopic level, whereas the latter approach tries to connect the observed macroscopic anisotropy of real porous materials with their microscopic geometry and structure. Permeability anisotropy caused by a finite size of a model system has not been regarded as an important factor so far, although some research techniques, e.g. numerical simulations, concentrate on artificially small porous systems. The reason for this lies in the fact that numerical flow simulations in complex porous geometries are extremely tedious and require extensive computer power and resources. A common strategy has been to perform calculations for just a few systems that are as large as possible Cancelliere et al. (1990); Başagaoglu et al. (2006). In contrast to this, here we solve the flow equations for hundreds or even thousands of different porous systems of small to medium sizes and then extrapolate the results to the limit of an infinitely large system. This method was already used in Matyka et al. (2008) to detect a small, systematic deviation of the flow tortuosity from several theoretical formulas, with an ad hoc interpretation of this phenomenon as a consequence of the finite-size anisotropy. Therefore, in this paper we present a systematic study of finite-size anisotropy in a two-dimensional model of statistically uniform porous media.
The structure of the paper is as follows. Section II specifies the model and the numerical techniques used. Main results are provided in Sec. III. Next, in Section IV we develop a simple theory to account for the asymptotic behavior of the angle between the external force and the volumetric fluid flux. Finally, the results are discussed in Sec. V.
In this study we use a model of freely overlapping squares Koponen et al. (1996, 1997); Matyka et al. (2008). In essence, this is a two-dimensional lattice model with a porous matrix modelled as a union of freely overlapping identical solid squares of size lattice units (l.u.) placed uniformly at random locations in a square lattice l.u. (). The squares are fixed in space but free to overlap, and their sides coincide with the underlying lattice. The remaining void space is filled with a fluid. A constant, external force is imposed on the fluid to model the gravity and we allow an angle () between the force and the system side to be arbitrary (note that in Koponen et al. (1996, 1997); Matyka et al. (2008) only the case was considered). Periodic boundary conditions are imposed in both directions to minimize finite-size effects. The porosity () is calculated as the ratio of unoccupied lattice nodes to the entire system volume (). The flow equations are solved in the creeping flow regime using the Lattice Boltzmann Model (LBM) Succi (2001) with a single relaxation time collision operator Bhatnagar et al. (1954) (see Matyka et al. (2008) for implementation details).
The model has three adjustable parameters: , , and . The first one corresponds directly to the macroscopic porosity. The value of affects the percolation threshold , which is a decreasing function of from (the standard site percolation threshold, ) Sahimi (1993) to (the continuous percolation threshold of aligned squares, ) Baker et al. (2002). As the model is solved using the LBM method without a numerical grid refinement Matyka et al. (2008), the minimum value of is 4 (this is the minimum length scale for the LBM method to resolve the macroscopic Navier-Stokes equations Succi (2001)). The value of controls the finite-size effects through the dimensionless ratio , which should be as small as possible to mimic an infinite system.
Anisotropy of fluid flow in the above-defined model will be investigated through Darcy’s law Bear (1972)
where is the volumetric fluid flux, is a symmetric tensor of the hydraulic conductivity, and is the unit vector in the direction of the gravitational field.
Iii Numerical results
iii.1 Tests on
The basic concept of transition from microscopic laws of hydrodynamics to macroscopic laws of transport in porous media is the representative elementary volume (REV), i.e. the smallest volume such that a measurement over it will yield a value representative of the whole Bear (1972). Darcy’s law (1) is, in principle, applicable only to systems that are larger than an REV, whereas significant anisotropy is expected in systems smaller than an REV. Hence, the primary question is whether or not Eq. (1) can be used to study anisotropy in small-size systems. To answer this we performed several simulations on , with its elements computed from for and .
First, the symmetry of was examined by quantifying the value of a dimensionless parameter given by
Furthermore, by choosing l.u., l.u., and , eighty different (i.e., fourty systems for each ) statistically uniform porous systems were constructed, for which 0.5% was found. This ensures that is symmetric within numerical errors in its elements. In a subsequent analysis we enforced -symmetry via replacing its off-diagonal elements ( and ) by their arithmetic mean, which ensures that is diagonizable.
Second, the tensorial properties of were examined by checking whether Eq. (1) can be used for an arbitrary . In particular, this equation predicts that if the mean flow direction () is aligned with the -axis, the angle between the external force () and the -axis should satisfy
This relation was examined for several systems, of which one is shown in Fig. 1, wherein, two streamline patterns for the same system ( l.u., l.u., ) with different -orientations are visualized. In the left panel, the external force is parallel to the (horizontal) -axis (), resulting in an angle of between the vector of the specific discharge () and the axis.
In the right panel, a force of the same magnitude makes an angle computed from (3); as expected, the angle between and the axis practically vanishes ().
iii.2 Tests on
A natural measure of anisotropy for a particular porous system is the angle between the vectors and . As this angle depends on the orientation of , following standard procedures in computer simulations, we fix . We verified that in this case the numerical value of (angle between and ) satisfies , which follows from symmetry arguments, and then calculated
In the above equation, denotes an average over different random porous systems. The results for l.u., , and several system lengths are shown in Fig. 2. For l.u. the data were fitted to
which yielded for and for . This suggests , i.e.
This relation does not hold for small ( in Fig. 2), for some realizations of such systems are likely to exhibit extreme anisotropy with so large that cannot be approximated by (for the angle between and can be as large as ).
Next we investigated statistical distribution of -values in different random systems with fixed , , and . In all cases this distribution closely resembled the normal distribution . Qualitative verification of this conjecture is presented in Fig. 3, which depicts the empirical cumulative distribution function (CDF) for two different system sizes (small symbols) together with the corresponding theoretical CDFs of the normal distribution,
The numerical data are in good agreement with (7). A quantitative comparison of the -distribution with was performed using the Kolmogorov-Smirnov test (at confidence level 95%). Out of all data points shown in Fig. 2 only the one corresponding to and did not pass the test, which in part is due to extremely large number of different samples (2000) used.
As mentioned before, the value of determines the percolation threshold for small porosities, and can be considered as a relevant parameter independent of and . For porosities much larger than , however, the connectedness and overlapping of individual randomly generated solid squares becomes irrelevant. In this case, using scaling arguments, one can expect that and are the only relevant parameters. Mathematically, this can be formulated as a similarity ansatz:
where is a similarity function. A direct numerical verification of this conjecture is difficult, as it requires averaging over many independent samples, which is rather a time-consuming job for large . Instead of this, we concentrated on a single parameter set with and that led to results demonstrated in Fig. 4
shown as cross symbols. These data were fitted to an ad-hoc formula
with three adjustable parameters , , and . The best-fit value of indicates that the approximation (8) can be safely used for .
Finally, we investigated the dependence of on porosity. One expects that should decrease from at the percolation threshold (a single, randomly oriented conducting channel) to at (completely permeable system). As shown in Fig. 5, our numerical results generally agreed with this picture. However, did not converge to its limiting value as . Instaed, it saturated at a positive value, which is independent of . Due to this rather unexpected result, we ensured that neither discretization errors nor large relaxation times affect the numerical data obtained for large porosities. We also verified that the system size used, , is sufficiently large for relation (6) to hold. This is clearly seen in the inset of Fig. 5, which depicts the product for and . The data for different collapsed in a broad range of . Porosities less than are in a vicinity of the percolation critical point, at which is expected to converge to as , and hence the product should diverge at . As the system size , it is possible that the the porosity range, for which scaling relation (6) does not hold, diminishes according to a power law. This behavior is a typical finite-size effect near a critical point Hunt (2005).
To explore the reason why does not tend smoothly to 0 as approaches 1, we inspected the streamlines in high-porosity systems exhibiting large anisotropy. An extreme example of such a system, generated with , is shown in Fig. 6a. At this high porosity, overlapping of individual obstacles is negligible, and the solid part of the system is made up of separate islands (that could correspond, for example, to a cross-section of a porous medium made of parallel fibers Koponen et al. (1998); Araújo et al. (2006)). Because the obstacles were placed uniformly and randomly in the whole system, their local concentration varies, and they form several larger groups of obstacles with relatively small distances between group members. Since fluid flux through a 2D channel is proportional to its width squared, most of the fluid flow takes place in relatively wide ‘channels’ between the groups. In other words, owing to no-slip boundary conditions on the obstacle surfaces, the fluid passes most easily in the inter-connected regions of low local obstacle concentration (high local porosity), whereas the regions of high local obstacle concentration (low local porosity)—even if occupied by separate obstacles—act effectively as large, almost impenetrable barriers. This solid-fluid ‘repulsion’ effect is not present in electric current flows (for the current intensity is proportional to the first power of a conductive channel width). For this reason, a high-porosity system which is highly anisotropic to fluid flow ( in Fig. 6a) exhibits a marginal anisotropy to electric current flow () as depicted in Fig. 6b.
iii.3 Tests on principal values
Mathematically, a porous system is anisotropic to flow if and only if at least two of the principal values of are diffrent. In the present case has two eigenvalues (principal permeabilities) and which can be ordered such that . Their ratio,
is equivalent to the ratio of the minimum to maximum permeabilities of a given porous system, and hence is a proper measure of its anisotropy Scholes et al. (2007). The more deviates from , the more anisotropic the system is.
We first verified that, as expected, the angle between the main principal axis and the -axis was distributed uniformly in the range (data not shown). Then the CDF of was determined for a particular case with , l.u., and l.u. As can be seen in Figure 7, the distribution of can be quite well fitted to the normal distribution . However, this is only an approximation, as in the present case for .
Iv Distribution of for large
Consider a porous system of size l.u. subject to an external force along the axis. Let denote the total displacement of a fluid particle as it passes the system between the opposite boundaries. While the -component of is a constant (equal to the system size ), the component (which we shall call lateral displacement and denote ) varies for different streamlines. If we calculate the average over all fluid particles, the angle between the volumetric fluid flux and the axis will satisfy
If is sufficiently small, this equation simplifies to
Let us consider a porous system of size l.u. and porosity . As shown in Fig. 8, it can also be regarded as two subsystems of size l.u. (labelled , ) or four subsystems of size l.u. (labelled ). Each of the small subsystems has its own permeability tensor , volumetric fluid flux , angle between the axis and , and mean lateral displacement with . Since the distribution of obstacles is uniform, porosities of each small subsystem is approximately equal to , and the mean lateral displacements can be considered as independent random variables drawn from the same distribution. Subsystems 1 and 2 form layer orthogonal to the external force. One may expect that the fluid streams passing through subsystems 1 and 2 are approximately the same in magnitude, and so the mean lateral displacement of the fluid, as it passes through layer , can be approximated by
Similarly, the mean lateral displacement in the layer can be approximated by . The mean lateral displacements of the fluid in layers and are practically independent of each other. This can be justified by an example of soil made of several horizontal and anisotropic layers—in this case the mean flow direction in a layer will depend only on the permeability tensor of this layer. This implies that the total lateral displacement of the fluid in the whole system () is approximately given by
where is calculated for the whole, system. Assuming that are independent random variables drawn from the same distribution with mean 0, one arrives at
which immediately leads to (6).
Equation (15) can be used iteratively to obtain
where the arguments of and (i.e. and ) indicate the system size. The right-hand side of this formula is an arithmetic mean of independent random variables with finite mean and variance, and—due to the central limit theorem—converges to normal distribution as . This explains why the distribution of for a sufficiently large system size can be approximated by a normal distribution (see Fig. 3).
The above can be readily extended to flows in an arbitrary space dimension . We skip the details and report only the final conclusions. First,
for sufficiently large . Second, the distribution of tends to the normal distribution as .
Equation (18) implies that anisotropy effects diminish with system size most quickly in three-dimensional (3D) systems (). Note, however, that the most important factor in computer simulations is the total number of lattice nodes (or volume) in the system. Using this quantity, equation (18) can be written as
irrespective of . Thus, anisotropy of the permeability tensor should be equally important (and difficult to account for) in computer simulations carried out in any space dimension.
It is important to verify Eq. (18) for space dimensions . While at the moment our software cannot be used for , the case can be tackled by studying a quasi one-dimensional system of size with fixed and . Note that in this case Eq. (18) should hold irrespective of whether the longer side of the system is parallel or perpendicular to the external force. The results, obtained for l.u., , , and ranging from to 800 l.u. are shown in Fig. 9 and confirm the validity of Eq. (18).
where is a function. This relation can be expected to hold in general only if and is sufficiently far away from the critical porosity . In a general case, is to be interpreted as a characteristic system length (such as the diameter of discs, in case the porous matrix is made of discs rather than squares), and depends on the system in question.
V Discussion and conclusions
Our results show that permeability anisotropy in statistically uniform porous systems of sizes typically used in computer simulations is a significant factor. The main parameter controlling this phenomena, especially at high porosities, is the ratio . For the model considered here, the asymptotic regime is observed for . In this regime the distribution of the angle between the external force and the volumetric fluid flux is very close to Gaussian, with the standard deviation diminishing as .
Although this conclusion is based on numerical results obtained for a particular model of a two-dimensional flow, it is expected to apply to a wide class of porous systems with randomly distributed identical solid matrices, such as squares, disks or spheres. This observation can be used to estimate the anisotropy-related statistical error in large-scale simulations, where often only one large system is considered for each parameter set Cancelliere et al. (1990). To this end it is enough to perform many independent simulations in small- and medium-size systems, verify that , and extrapolate to the required value of . Next, assuming that the distribution of is normal, one obtains the complete information about the error related to the anisotropy of the permeability tensor.
Magnitude of permeability anisotropy could serve as a good indicator of how the size of a model system compares with that of a REV. We found that even for the angle between the external force and the volumetric fluid flux can be as large as , and the permeability can vary with the orientation of the external force by a factor of 2. The value below which the anisotropy effects are small enough to be practically negligible is , as in this case , i.e. with probability . This enables to estimate the size of a REV in the model considered here as l.u.
It is interesting to note that most of the simulations carried out so far for 2D systems do not meet the criterion of , mainly because they used models with large . In previous studies on two-dimensional flows in various statistically uniform porous media, many researchers used -values ranging from J. S. Andrade, Jr. et al. (1995), through Verberg and Ladd (1999), 0.04 Kostek et al. (1992); Martys and Garboczi (1992), 0.05 Matyka et al. (2008); Yiotis et al. (2007) to 0.1 Koponen et al. (1996, 1997); Nabovati and Sousa (2007), usually assuming their systems to be isotropic. In view of our present findings, validity of this assumption in some of these cases is questionable and requires verification. Generally, one should expect that the threshold value of below which the permeability anisotropy is negligible is not universal, but depends on the geometry and structure of the porous medium, especially on its porosity and space dimensionality.
Anisotropy is a phenomenon independent of the boundary conditions. Periodic boundary conditions used in this paper facilitate measurement of the permeability tensor and reduce finite-size (boundary) effects. Other boundary conditions could mask, but would not eliminate anisotropy effects. For example, using solid walls along the fluid flow would fix the orientation of the fluid stream, however, the system would respond to such boundaries with an internal pressure gradient Bear (1972), which would change (and complicate measurement of) the orientation of the effective force acting on the fluid.
Finite-size permeability anisotropy in three-dimensional small porous systems remains an open problem. Typical system sizes used in numerical 3D simulations are l.u. The ratio is thus much larger in 3D than in 2D simulations and ranges from Knackstedt and Zhang (1994); Zhang and Knackstedt (1995), through Yiotis et al. (2007), Succi et al. (1989); Cancelliere et al. (1990), to 0.33 Verberg and Ladd (1999). The magnitude of permeability anisotropy is usually neglected. One exception is the paper by Verberg and Ladd Verberg and Ladd (1999), who published the off-diagonal elements of the permeability tensor. Their data for a single configuration of randomly distributed spheres suggests that is a decreasing function of the porosity and varies from for to for . This is in agreement with our present findings for a 2D system and indicates that permeability anisotropy is especially important close to the percolation threshold.
Acknowledgements.We gratefully acknowledge support from UWr grant Nr 2944/W/IFT/08 (ZK, MM), Sparda Bank Münster eG (MM), and Max Planck Institute for Marine Microbiology in Bremen (MM).
- Martys and Garboczi (1992) N. Martys and E. J. Garboczi, Phys. Rev. B 46, 6080 (1992).
- Koponen et al. (1997) A. Koponen, M. Kataja, and J. Timonen, Phys. Rev. E 56, 3319 (1997).
- Nabovati and Sousa (2007) A. Nabovati and A. C. M. Sousa, J. Eng. Sci. Techn. 2, 226 (2007).
- Verberg and Ladd (1999) R. Verberg and A. J. C. Ladd, Phys. Rev. E 60, 3366 (1999).
- Kostek et al. (1992) S. Kostek, L. M. Schwartz, and D. L. Johnson, Phys. Rev. B 45, 186 (1992).
- Matyka et al. (2008) M. Matyka, A. Khalili, and Z. Koza, Phys. Rev. E 78, 026306 (2008).
- Başagaoglu et al. (2006) H. Başagaoglu, P. Meakin, S. Succi, and J. Welhan, Europhys. Lett. 73, 858 (2006).
- Zhang and Knackstedt (1995) X. Zhang and M. A. Knackstedt, Geophys. Res. Lett. 22, 2333 (1995).
- Zerai et al. (2005) B. Zerai et al., Transp. Porous Med. 60, 159 (2005).
- Bear (1972) J. Bear, Dynamics of Fluids in Porous Media (Elsevier, New York, 1972).
- Scholes et al. (2007) O. N. Scholes, S. A. Clayton, A. F. A. Hoadley, and C. Tiu, Transport in Porous Media 68, 365 (2007).
- Boudreau and Meysman (2006) B. P. Boudreau and F. J. Meysman, Geology 34, 693 (2006).
- Alam et al. (2006) P. Alam, T. Byholm, and M. Toivakka, Nord. Pulp Pap. Res. J. 21, 670 (2006).
- Ohkubo (2008) T. Ohkubo, Transp. Porous Med. 72, 339 (2008).
- Cancelliere et al. (1990) A. Cancelliere, C. Chang, E. Foti, D. Rothman, and S. Succi, Phys. Fluids A 2, 2085 (1990).
- Koponen et al. (1996) A. Koponen, M. Kataja, and J. Timonen, Phys. Rev. E 54, 406 (1996).
- Succi (2001) S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond (Clarendon Press, New York, 2001).
- Bhatnagar et al. (1954) P. L. Bhatnagar, E. P. Gross, and M. Krook, Phys. Rev. 94, 511 (1954).
- Sahimi (1993) M. Sahimi, Rev. Mod. Phys. 65, 1393 (1993).
- Baker et al. (2002) D. R. Baker, G. Paul, S. Sreenivasan, and H. E. Stanley, Phys. Rev. E 66, 046136 (2002).
- Hunt (2005) A. G. Hunt, Percolation Theory for Flow in Porous Media (Springer, Berlin, 2005).
- Koponen et al. (1998) A. Koponen, D. Kandhai, E. Hellén, M. Alava, A. Hoekstra, M. Kataja, K. Niskanen, P. Sloot, and J. Timonen, Phys. Rev. Lett. 80, 716 (1998).
- Araújo et al. (2006) A. D. Araújo, W. B. Bastos, J. S. Andrade, Jr., and H. J. Herrmann, Phys. Rev. E 74, 010401(R) (2006).
- J. S. Andrade, Jr. et al. (1995) J. S. Andrade, Jr., D. A. Street, T. Shinohara, Y. Shibusa, and Y. Arai, Phys. Rev. E 51, 5725 (1995).
- Yiotis et al. (2007) A. G. Yiotis, J. Psihogios, M. E. Kainourgiakis, A. Papaioannou, and A. K. Stubos, Colloid Surf. A–Physicochem. Eng. Asp. 300, 35 (2007).
- Knackstedt and Zhang (1994) M. A. Knackstedt and X. Zhang, Phys. Rev. E 50, 2134 (1994).
- Succi et al. (1989) S. Succi, E. Foti, and F. Higuera, Europhys. Lett. 10, 433 (1989).