Investigation of pore-scale random porous media using lattice boltzmann method

1.

Introduction

Fluid flow through permeable materials is one of the topics of interest in the geosciences, petroleum industry, molding industry and etc. In all above mentioned industries we are faced with flows with relatively low Reynolds numbers where Darcy’s law is applicable. Permeability is the most important property which characterizes a porous medium. It is a measure of the frictional resistance of the material to fluid flow or, equivalently, the drag force of the fluid on the material. Hence, it needs to be specified prior to any macro scale modelling.

Two major groups of materials can be assumed as an appropriate representation of porous media. The first one is a fibrous medium. A fibrous medium is composed of a number of cylinders which are called fibers. The axes of fibers can be parallel to each other (1D), placed on planes parallel to each other with random orientations in each plane (2D), or randomly oriented in space (3D). The second group is granular media, and as a representative, an array of spheres which are considered as solids in a porous medium.

Many experimental works have been carried out for granular media to model the Darcy and Darcy–Forchheimer drags: the Blake–Kozeny [1] equation, the Ergun equation [2] and any modifications thereof [3, 4]. Sangani and Acrivos [5] studied the drag in periodic arrays of spherical particles; Liu et al [6] proposed a semi-empirical formula for the pressure drop which incorporates the tortuosity, the curvature ratio and the variation of the pore cross-sectional area. For 3D flows, Larson and Higdon [7] performed calculations for the Stokes flow through a lattice of spheres. Kim and Russel [8] derived an expression for drag of random arrays of spherical particles. Nakayama et al [9] studied the flows through a lattice of cubes to estimate the contributions of both Darcy and Forchheimer drags to the macroscopic pressure drop.

Treating the flow through porous medium by the method of the conventional Navier–Stokes codes is often faced with extensive computational time, poor convergence and/or numerical instability that stem from the narrowness of the flow passages. In recent years, the Lattice Boltzmann Method (LBM), as a mesoscale approach, has emerged as an alternative tool to investigate flow in complex geometries. In general, the Lattice Boltzmann Method (introduced in the next section) is easier to be implemented than conventional computational fluid dynamic techniques, is highly compatible with parallelization, and can deal with arbitrary complex flow geometries without heavy penalties. These attributes, combined with its flexibility, make LBM suitable for the investigation of a wide range of flow problems ranging from flows through complex geometries (e.g. permeating groundwater flow, melt segregation) to flows in which multiphase interactions are of interest (e.g. particle laden flows, bubble suspensions). Rothman [10] and Chen et al [11] used the lattice gas automata to study the microscopic behavior occurring at the pore scale and obtained volume-averaged parameters from a microscopic point of view. Succi et al [12] and Cancelliere et al [13] employed the LBM to extract the permeability of a randomly distributed 3D porous medium. Using lattice-Boltzmann approach, Maier et al. [14] studied flow and transport properties in packed columns of spheres. Inamuro et al [15] applied the LBM to examine flows through a 3D porous structure, which was composed of nine identical spheres in a rectangular domain, for high and low Reynolds numbers. Their results were in good agreement with the Erguns correlation. However, they covered only one case of porosity and structure. Manwart et al [16] carried out a comparative study on LBM and a conventional FDM for flows through a straight rectangular channel and a cubic array of spheres with a porosity of 0.15. Each method was evaluated using the exact solutions of the Stokes equations. Next, both algorithms were employed to estimate the permeability of the three-dimensional sandstones. Van der Hoef et al. [17] investigated the drag in an arrangement of spherical particles with binary size distributions. The effect of particle shape on permeability has been studied numerically by Coelho et al. [18], Stewart et al. [19], and Garcia et al. [20]. All the above mentioned studies are restricted to the continuum flow regime even though the LBM can be applied to the slip flow regimes in micro-scale, and despite the seemingly large volume of experimental or computational literature, a detailed study on various porous structures, shapes of composed materials and porosities remains scanty. Jeong et al. [21] studied the macroscopic porous medium of various structures by calculating flows through 2D and 3D porous structures. Moreover, they investigated the effect of Knudsen number, Kn, on macroscopic porous-medium properties.

In this study, pore-scale transport properties in the porous medium of various structures are investigated for flows through 2D and 3D porous structures, including 2D unit-cell with randomly distributed rectangular obstacles and 3D porous materials composed of randomly distributed spheres. The effect of geometrical parameters of media on permeability and tortuosity is also studied. The results are validated by various analytical, experimental, and numerical data from literature. Simulations are performed using LBFlow package source codes. Code is available by request at http://www.dur.ac.uk/ed.llewellin/lbflow/. It has been modified by authors.

2. Numerical Method

2.1 D2Q9 and D3Q15 Lattice Boltzmann Scheme

The LBGK[a] [22] lattice Boltzmann scheme is governed by

with  the particle distribution in the direction  ,  the position,  the discrete velocity in the direction  , the time step,  the dimensionless relaxation time,  the fluid viscosity, and  the lattice sound speed, equal to .

In the present study, the two-dimensional 9-velocity (D2Q9) and the three-dimensional 15-velocity (D3Q15) lattice Boltzmann models are used [23]. Schematic diagrams of D2Q9 and D3Q15 models are shown in Figure 1.

In these models, the discrete velocities  are given by

The equilibrium distribution function, , is expressed by [24],

is the weighting factor of the system, is the density, and  is the macroscopic velocity. The values of  are given as [25]

In LBM, equation (1) is divided into two steps: collision and streaming.

The collision step is:

and streaming step is :

is the post-collision state of the distribution function . After the simulation of the single variable , the macroscopic flow properties, such as the density, pressure, and momentum, are extracted as follows:

2.2 LBM Dimensional Consideration

Physical units of length, time and mass are distinct from the simulation units used internally by the lattice Boltzmann algorithm. For a simulation of a practical usage, we require a unit conversion between simulation and physical units [26]:

Equations (10) are usually called mapping relationships and are used for conversion of physical and lattice units to each other. , , , , , and  are lattice length, lattice time, lattice mass, and mapping parameters for length, time, and mass, respectively. The mapping allows us to convert, for instance, velocities from lattice spacing per time-step in the simulation, to meters per second in physical system. In this study, mapping parameters are entered to the simulator, which will be discussed later.

3. Boundary Conditions and Geometry Generation

3.1 Boundary Conditions

One of the advantages of LBM is its easy boundary condition implementation on complex geometries including complex walls [24]. The following three boundary conditions are used for different conditions in the present study:

3.1.1 Halfway Bounce-back BCs :

The halfway bounce-back boundary condition is applied on a solid surface [23]. Figure 2 illustrates the halfway bounce-back scheme [23].

Fundamental steps of LBM, namely collision and streaming^[b] are completely local; hence they are done in a fashion regardless of the type of the geometry. That is, there is no difference between the ways that collision is done in a random or a regular medium. There are two kinds of nodes in a porous medium: solids and fluids.

Streaming and collision are both applied only for fluid nodes. Now if it is assumed that collision step is done for all fluid nodes including those next to solid boundary nodes, a problems is just encountered when it’s time for the streaming step to be done; It is mentioned that collision and streaming are only valid for fluid nodes, hence no collision/streaming rules exist for solid boundary nodes. If so, what happens for a fluid node, just next to boundary solid nodes? In other words, in figure 2 bellow, because solid nodes are not taken part into collision/streaming, upcoming f₂, f₅ and f₆ of the adjacent-to-boundary fluid nodes are unknown. Consequently, the calculation becomes unsettled.

Here’s when boundary conditions are employed. Boundary Conditions are special kinds of collision steps which are set to solve the problem of unknown streaming data of adjacent-to-boundary fluid nodes. In this example the bounce-back boundary condition is employed as : f₅₌ f_7, f₂₌ f₄, and f₆₌ f₈ to make unknown distribution functions for adjacent-to-boundary nodes specified.

The solid boundary is located at the halfway between fluid node and solid isolated node. , and  are known from streaming process, but , and are unknown distribution functions as

they come from isolated solid nodes outside the computational domain.

3.1.2 Pressure BCs (Dirischlet):

Pressure boundary conditions were applied along the x direction, which was supposed to be the major flow direction. This kind of boundary condition constrains the density[c] at the boundary. Two densities were specified at inlet and outlet of the domain along with the major flow direction. To make the point clear, consider an entrance boundary node as illustrated in figure 3. In this node, there are known and unknown distribution functions, and unknown normal velocity, :

Note that velocity tangent to the boundary, u_yis assumed to be zero. According to above figure, there are four unknowns. Hence, a system of four equations is needed.

Three equations are written according to the moment equations, Eqs. 9:

For the fourth equation, Zou and He [44] suggested that non-equilibrium components of  and  are assumed equal:

Equations 11 and 12 together are a system of four equations with four unknowns. By solving this system, four unknowns, , and , are determined. The same procedure is performed for all inlet and outlet boundary nodes.

3.1.3 Periodic BCs :

Periodic boundary conditions are applied for the transverse directions[d] in which the system becomes closed if one edge is attached to the opposite edge. As an example, in figure 4 for inlet and outlet boundaries:

where and  denote inlet and outlet, respectively.

3.2 Geometry Generation

Two-dimensional domains with a mesh size of 100×100 lu² and three-dimensional domains with a 100×100×100 lu² mesh size were selected with periodic boundary conditions in all three spatial directions. To induce the flow, pressure difference was applied between the inlet and the outlet boundaries [21].

In order to reconstruct the required geometries, instead of tomography of real porous media, an alternative method, namely computed tomography was employed. As the input geometry for 3D simulations, some 3D virtual media, for example

randomly distributed spheres in a unit cell, were generated using C++ coding language. Next, each volume was parsed into a sequential stack of 2D bmp image files, each file as a representation of a slice of the volumetric medium. Then the resulting stack of images was imported into an image processing program to be converted to an identifiable format for the simulator. Another geometry generation method was to create some ASCII text files, including only “0”s and “1”s as fluid and solid nodes, respectively. Generated geometries were put to the simulations. Figures 5 and 6 illustrate binary and ASCII geometry files. For 2D medium which has rectangular obstacles, the aspect ratio is the ratio of the height to the length of obstacles. For example in the medium depicted in figure 5, right, the aspect ratio of obstacles is 0.5.

Details of geometry generation are included in appendix.

For engineering purposes, the permeability, K, is generally non-dimensionalized by dividing it by the square of the characteristic length scale. For a pore-scale simulation, it is customary to use the hydraulic diameter of obstacles, D_h, as the characteristic length.

where  and  denote number of dimensions, area, and perimeter of obstacles, respectively.

Now the dimensionless permeability is defined as  :

4. Results and Discussion

4.1 Two-dimensional Simulations

After putting the input geometry file into the simulator, some parameters have to be entered as well, including the pressure difference in  direction, the kinematic viscosity , the lattice size, the relaxation parameter,  , and mapping parameters[e]. Boundary conditions are implemented in the main CPP file. The simulation is started and the mean superficial velocity (i.e. Darcy velocity) is calculated as the main output by applying the Darcy’s law:

It should be noted that when the assumption of isotropic medium is not possible, the general form of Darcy's law should be considered, which is given in equation 17. This general form is reduced to equation 16 for assumed homogenous and isotropic porous media.

As mentioned before, the pressure difference is applied in x direction, and the medium is assumed to be isotropic. Hence, equations 16 and 17 are reduced to equation 18:

By the use of this equation, the permeability is extracted. Broadly speaking, for engineering purposes, the permeability K is non-dimensionalized by using equation 15. It is worth noting that the Darcy’s law is valid while the Reynolds number value is very low[f]. In order to keep the Reynolds number low, the value of applied pressure gradient has to be absolutely low.

Figure 7 demonstrates contour of velocity magnitude in 2D random geometries for creeping flow regime as investigated in the current study. Contours are generated using LB2D_PRIME open-source simulator.

4.1.1 Permeability

Results by LBM are evaluated against the results of Carman-Kozney [27] and modified Carman-Kozney correlations presented by Koponen et al. [28]. Note that the dimensionless permeability is very sensitive to changes in the porosity, and the graph is plotted with the ordinate in logarithmic scale. Koponen et al. [28] used effective porosity, , based on neglecting dead-end pores. Then, they proposed the following correlation, relating the effective porosity to real porosity of the medium,

They also reported correlations for hydraulic tortuosity and specific surface, S, dependent to porosity [28, 29]. Substituting their correlations for tortuosity and specific surface in their permeability correlation leads to their dimensionless permeability correlation.

If one uses real porosity, the Carman-Kozney correlation is concluded, which is proposed in [27]. Results are shown in figure 8. Scattered data are the results of present work for 73 random 2D samples with different porosities and obstacle aspect ratios. In addition, a curve with fitting parameters with the minimum error is fitted over scattered data, leading to dimensionless correlation for permeability.

with fitting parameters equal to 0.001, -0.268, and -19.95 . The “correlation coefficient”, , and the “ coefficient of determination”, , of this correlation have excellent values of 0.997 and 0.994, respectively. The overall error of this fitted correlation is 1.63%. Standard deviation,  , for this case is 2.316292. This correlation is plotted as red dashed-line in figure 8. From figure 8, it can be inferred that the results of this study are in good agreement with the results of Koponen et al. [28] for porosities higher than 0.65. However, while the porosity is decreased to less than 0.65, some deviations are emerged, that is to say, the LBM overestimates the permeability in comparison with that of Koponen et al. [28]. It is likely because of the increase of dead-end pores with the decrease of porosity; hence the effective porosity is decreased. Subsequently, a rise in the difference between real and effective porosities occurs that leads to a sharp fall in graph of Koponen et al. [28] correlation with porosities less than 0.65.

The results of dimensionless permeability according to different obstacle aspect ratios are depicted in figure 9. Three different obstacle aspect ratios are employed: 0.5, 1, and 2. From this figure, for aspect ratios equal to 0.5 and 1 in porosity range of   there is some mild data variation. However, a smooth increase rate in comparison with aspect ratio equal to 2. Thus, it can be concluded that a choice of 0.5 and 1 values for aspect ratio may lead to more accurate results.

4.1.2 Tortuosity

The hydraulic tortuosity is expressed as [28]:

where L_eff is the real length of the flow through interconnected pores of the porous medium, and  is the minimum available length, as if the medium is not porous. A general discussion in Ref. [30] leads to the following form of tortuosity:

where is the average magnitude of the intrinsic velocity over the entire volume and is the volumetric average of its component along the macroscopic flow direction.

Here, the fact that the LBM method uses a regular mesh allows to approximate Eq. 22 with:

where  runs through all lattice nodes. Note that this simple formula can be used not only in numerical studies, but also for the data obtained experimentally.

Figure 10 depicts the effect of the obstacle aspect ratios on the predicted tortuosity. As it is observed, the flow tortuosity increases with increasing obstacle aspect ratios. The effect of obstacle aspect ratios is more pronounced at lower porosities, however at high porosities it is practically negligible. Another fact extracted from graphs in figure 10 is that the lower is the porosity, the greater would be the tortuosity values. Physical interpretation of this fact may be as follows: while the porosity is decreased, the total volume of dead-end pores is increased, making the fluid to do more trial and errors to find interconnected pores, leading to a longer fluid’s travelling distance, and ending to a more tortuosity value based on Eq. 21.

4.1.3 Required Time-steps

From figure 10 one may conclude that tortuosity is not only a function of porosity, but also a function of obstacle aspect ratios.

In order to define a criterion for convergence of the solutions, the Convergence Number is defined as follows:

Whenever the deviation of results of 50 sequential time-steps is less than Criteria Number, the simulation will be terminated.

The choice of Convergence Number is done by trial and error. Table 1 depicts eclectic Convergence Numbers for diverse porosity ranges.

The effect of porosity on required time-steps for convergence of the results for different obstacle aspect ratios is reported in figure11. From graphs, it can be inferred that:

i. Convergence speed is increased through the middle of the porosity range.

ii. With an increase in porosity, number of required time-steps for convergence is increased. Since the increase of the porosity is equal to the rise of number of fluid nodes, and this consequently leads to more collisions and propagations, which are two major LBM mechanisms, more local macroscopic velocities are to be calculated by eq. 9. That is to say, the total number of time-steps is increased.

iii. On the other part, figures clearly show that a rise of the number of time-steps is also occurred with a reduction in the porosity. It is because of the fact that at low porosities the volume of dead-end pores is increased, which increases the total fluid travel distance and consequently more processing time.

Hence, the time-step is increased along two end limits of the porosity band, and a sharper increase of time-step at high porosities is clearly observed in all three graphs in figure 11.

4.2 Three-dimensional Simulations

4.2.1 Permeability

Dimensionless values of extracted permeability in some samples of 100×100×100 lu³sphere-packing, versus porosity are shown in figure 12. Results are validated using three prevalent correlations:

               i.            Kozney-Carman [27],

              ii.            Rumpf-Gupte [31], and

            iii.            Koponen et al. [28].

Scattered data are the results of calculation of permeability over 29 3D random granular media. Similar to the 2D case, a correlation is proposed based on the least-squares method by fitting the parameters with minimum fitting errors:

with fitting parameters equal to 0.001, -0.268, and -19.95. The “correlation coefficient”,  and the “coefficient of determination”, of this correlation have excellent values of 0.995 and 0.99, respectively. The overall error of this fitted correlation is 4.28%. And the Standard deviation,  for this case is 2.55136. The red dash-line depicts correlation 24.

Some good agreements between the data can be seen for moderate to high values of porosity. However, there are some deviations in upper and lower limits of porosity.

To investigate the effect of spheres’ diameter on results, three distinct scattered graphs corresponding to three different sphere diameters are depicted in figure 13. It can be seen that when the lower is the sphere diameter, the higher would be the calculated permeability at a constant porosity.

4.2.2 Tortuosity

Figure 14 illustrates the results of 3D tortuosity which are validated using two experimental, and one analytical correlations. Various phenomenological expressions have been proposed to describe the tortuosity as a function of the porosity. Among them, the logarithmic equation is valid for media with non-porous, non-overlapping particles [32].

There is another experimental power-law correlation which is prevalent for granular media, again with non-overlapping spheres [33-36] with a coefficient of , being dependent on type of sphere-packing. According to reports of Refs. [36-41], the value of  for non-overlapping condition would be 0.4

The third correlation is an analytical one, namely the Maxwell equation[g] [42]. A comparison between the results is depicted in figure 14. The figure clearly shows that the results of the present study agrees well with the results of Maxwell’s analytical correlation. However, there are some drastic deviations from the experimental results, which are based on the assumption of overlapping spheres in this study and non-overlapping assumption in the experiments. Hence, the present study data, based on overlapping spheres, are fitted to experimental and analytical non-overlapping formulas, with fitted parameters listed in table 2.

5. Conclusions

The flow behavior, permeability, and tortuosity in two and three-dimensional random porous media were investigated using Lattice Boltzmann Method with D2Q9 and D3Q15 lattice arrangements. Unit-cells with randomly distributed rectangular obstacles, and random sphere-packing in cubic unit-cells were employed as two and three-dimensional media, respectively. Almost all models, including LBM models and models available in the literature behave the same at low porosities: no significant difference in permeabilities is seen at low porosities. However, sudden sharp rises at high porosities are observed. Thus, the graph was plotted with the ordinate in logarithmic scale. For two-dimensional cases, some deviations from the permeability results of Koponen et al. [8] at the low end of porosity were seen. It occurred because of the increase in the difference between effective and real porosities at low porosities.

The effects of obstacle aspect ratios and spheres’ diameter on permeability were also studied for two and three-dimensional cases. It is observed that the when the lower are the obstacle aspect ratios and diameters, the higher would be the calculated permeability at a constant porosity, and the more uniform would be the variation of the permeability. Based on the scattered results, two permeability  correlations  for  both  two  and  three- permeability correlations for both two and three-dimensional cases with the least fitting errors were proposed. Tortuosity was also studied for both cases and the effect of porosity and aspect ratio was depicted. It is concluded that tortuosity is a function of not only

porosity, but also obstacle aspect ratios. Finally, it is worth mentioning the effect of porosity on the required time-steps for solution to converged.

Appendix. Geometry generation in details

Two types of geometries were generated in this paper: Binary and ASCII.

1. Details of binary geometry generation:

These type of geometries simulate porous media by assuming obstacles as white color (grayscale=255) and pores as black color (grayscale=0).

2D binary geometries are simply created by Microsoft Paint application, drawing white rectangles on a black background as shown in figure A1:

Then the porosity of just generated 2D binary geometries is measured with ImageJ, an open-source image processing software[h].

For 3D binary cases, a simple C++ program was employed to generate random spheres in a cube as a unit  cell  domain.  The  code  simply  chooses  random geometrical positions as centers of spheres in the defined domain[i], then the C++ sphere generation algorithm is employed to create diversely placed spheres as random obstacles.

The output of C++ is a “.raw” format file which can be broken up to a sequence of “.bmp” files with ImageJ software, thus its porosity can be calculated. The resulted “.raw file” can be entered to the LBM simulator, that is to say LB3D-Prime simulator.

2. Details of ASCII geometry generation:

Another approach of geometry generation is to create geometries in a simple ASCII text file, with 0s and 1s as representations of pores and obstacles, respectively. Granted the fact that LBFlow, another simulator, uses ASCII geometries as input files, creating the ASCII is an indispensable step. Thus, after creating binary geometries with above approaches, a Netpbm script in Linux is written to convert the files into ASCII files: by putting all relevant “.bmp” files[j] into one directory, and then

by typing this script on Linux command line:
cat <(j=0; for i in fence*.bmp;

do j=$((j+1));

done;

echo "100 100 $j";

for i in fence*.bmp;

  do bmptopnm $i | pnmtoplainpnm |

  sed ':a;N;$!ba;s|P2\n100 100\n255\n||; s|255|1|g;

  s| ||g; s|\n||g';

done) > mask.dat

The script above is modified for a  domain. The expression “for i in fence*.bmp” asks the Linux compiler to place all “.bmp” files in a ASCII file named “mask.dat”. First, it is required that all “.bmp” files are converted to “.pnm” files, and then converted to a set of 0/1 binary characters, as is achieved through lines 6-8.

The resulting “.mask” file, figure 6, can now be entered into the numerical simulator.

References

[1].    R.B. Bird, W.E. Stewart, and E.N. Lightfoot, Transport Phenomena, Wiley, New York, P. 199 (1960).

[2].    Ergun S 1952 Fluid flow through packed columns Chemical Engineering Progress, 48, 89-94, 1952.

[3].    Macdonald I F, EL-Sayed M S, Mow K and Dullien F A L 1979 Flow through porous media—the Ergun equation revisited Industrial & Engineering Chemistry Fundamentals, 18 199-208, 1979.

[4].    R. M. Fand, B.Y.K Kim, A.C.C. Lam, and R.T. Phan, Resistance to the flow of fluids through simple and complex porous media whose matrices are composed of randomly packed spheres, Journal of Fluids Engineering, 109, 268, 1987.

[5].    A. S. Sangani, A. Acrivos, Slow flow through a periodic array of spheres, International Journal of Multiphase Flow, 8(4), 343-360, 1982.

[6].    S. Liu, A. Afacan, and J. Masliyah, Steady incompressible laminar flow in porous media, Chemical Engineering Science, 49, 3565, 1994.

[7].    R. E. Larson, and J.J.L. Higdon, A periodic grain consolidation model of porous media, Physics of Fluids A, 1, 38-46, 1989.

[8].    S. Kim, W.B. Russel, Modeling of porous media by renormalization of the Stokes equations, Journal of Fluids Mechanics, 154, 269-286, 1985.

[9].    A. Nakayama, F. Kuwahara, Y Kawamura, and H. Koyama, Three-dimensional numerical simulation of flow through a microscopic porous structure Proc. ASME/JSME Thermal Engineering Conference, 3, 313, (1995).

[10]. D.H. Rothman, Cellular-automaton fluids: a model for flow in porous media, Geophysics, 54, 509, 1988.

[11]. S. Chen, K. Diemer, G.D. Doolen, K. Eggert, C. Fu, S. Gutman, and B.J. Travis, Lattice gas automata for flow through porous media, Physica D, 47, 72-84, 1991.

[12]. S. Succi, E Foti, and F. Higuera, Three-dimensional flows in complex geometries with the lattice Boltzmann method Europhysics Letter, 10, 433-438, 1989.

[13]. A. Cancelliere, C. Chang, E. Foti, D.H. Rothman, and S. Succi, The permeability of a random medium: comparison of simulation with theory, Physics of Fluids, A 2, 2085-2088, 1990.

[14]. R.S. Maier, D.M. Kroll, Y.E. Kutovsky, H.T. Davis, and R.S. Bernard, simulation of flow-through bead packs using the lattice boltzmann method, Physics of Fluids, 10 (1), 60-74, 1998.

[15]. T. Inamuro, M. Yoshino, and F. Ogino, Lattice Boltzmann simulation of flows in a three-dimensional porous structure, International Journal for Numerical Methods in Fluids, 29, 737-748, 1999.

[16]. C. Manwart, U. Aaltosalmi, A. Koponen,  R. Hilfer, and J. Timonen, Lattice-Boltzmann and finite-difference simulations for the permeability for three-dimensional porous media, Physical Review E, 66, 016702, 2002.

[17]. M.A. van der Hoef, R. Beetstra, and J.A.M. Kuipers, Lattice-Boltzmann simulations of low-Reynolds-number flow past mono- and bidisperse arrays of spheres: results for the permeability and drag force, Journal of Fluids Mechanic, 528, 233-254, 2005.

[18]. D. Coelho, J.F. Thovert, and P.M. Adler, Geometrical and transport properties of random packings of spheres and aspherical particles, Physical Review, E 55, 1959-1978, 1997.

[19]. M.L. Stewart, A.L. Ward, and D.R. Rector, A study of pore geometry effects on anisotropy in hydraulic permeability using the lattice-Boltzmann method, Advanced Water Resources, 29, 1328-1340, 2006.

[20]. X. Garcia, L.T. Akanji, M.J. Blunt, S.K. Matthai, and J.P. Latham, Numerical study of the effects of particle shape and polydispersity on permeability , Physical Review E 80(2), 222, 145–197, 2009.

[21]. N. Jeon, D.H., Choi, and C.L. Lin, Prediction of Darcy-Forchheimer drag for micro-porous structures of complex geometry using the lattice Boltzmann method, Journal of Micromechanics and Microengineering, 16, 2240–2250, 2006a.

[22]. P.L. Bhatnagar, E.P. Gross, and M. Krook, A model for collisional processes in gases I: small amplitude processes in charged and in nutral one-component systems, Physical Review, 94, 511-525, 1954.

[23]. A.A. Mohammad, Lattice Boltzmann Method, Fundamentals and Engineering Application With Computer Codes, Springer, London, (2011).

[24]. S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, Oxford, UK, pp288, (2001).

[25]. Y.H. Qian, D. d’Humières, and P. Lallemand, Lattice BGK models for Navier–Stokes equation. Europhysics Letter, 17, 479–484, 1992.

[26]. E.W. Llewellin, LBflow: An extensible lattice Boltzmann framework for the simulation of geophysical flows. Part I: theory and implementation”. Computers & Geosciences, 36, 115–122, 2010.

[27]. P.C. Carman, Fluid flow through granular beds, Transactions of the Institution of Chemical Engineers, 15, 150–166, 1937.

[28]. A. Koponen, M. Kataja, and J. Timonen, Permeability and effective porosity of  porous media, Physical Review E, 56 (3), 3319- 3325, 1997.

[29]. A. Koponen, M. Kataja, and J. Timonen, Tortuous flow in porous media, Physical Review E, 54 (1), 406-410, 1996.

[30]. A. Duda, Z. Koza, and M. Matyka, Hydraulic tortuosity in arbitrary porous media flow, Physical Review E 84, 036319, 2011.

[31]. H. Rumpf, A.R. Gupte, The influence of porosity and grain size distribution on the permeability equation of porous flow, Chemie Ingenieur Technik (Weinheim), 43 (6) , 367-375, 1975.

[32]. E. Mauret, M. Renaud, Transport phenomena in multi-particle systems—I. Limits of applicability of capillary model in high voidage beds-application to fixed beds of fibers and fluidized beds of spheres, Chemical Engineering  Science, 52, 1807-1817, 1997.

[33]. F.G. Ho, W. Strieder, A variational calculation of the effective surface diffusion coefficient and tortuosity, Chemical Engineering  Science, 36, 253–258, 1981.

[34]. H. Pape, L. Riepe, and J.R. Schopper, Interlayer conductivity of rocks-a fractal model of interface irregularities for calculating interlayer conductivity of natural porous mineral systems, Colloids and Surfaces, 27, 97–122, 1987.

[35]. M.R. Riley, F.J. Muzzio, H.M. Buettner, and S.C. Reyes, A simple correlation for predicting effective diffusivities in immobilized cell systems, Biotechnology and Bioengineering, 49, 223–227, 1996.

[36]. M. Mota, J.A. Teixeira, and A. Yelshin, Binary spherical particle mixed beds: porosity and permeability relationship measurement, Transactions of the  filtration Society, 1, 102–106, 2001.

[37]. K. Klusáček, P. Schneider, Effect of size and shape of catalyst microparticles on pellet pore structure and effectiveness, Chemical Engineering  Science, 36, 523–527, 1981.

[38]. R.J. Millington, J.P Quirk, Permeability of porous solids, Transactions of the Faraday Society, 57, 1200–1207, 1961.

[39]. T.C. Zhang, P.L. Bishop, Evaluation of tortuosity factors and effective diffusivities in biofilms, Water Research, 28, 2279–2287, 1994.

[40]. A. Yelshin, M. Mota, and J. Teixeira, Proceedings of the International Conference on Filtech Europa-97, Dusseldorf, Filtration Society, Horsham, UK, October 14–16, 327–334 (1997).

[41]. M. Mota, J.A. Teixeira, and A. Yelshin, Ed. Feyo de Azevedo, E. Ferreira, K. Luben, O Osseweijer (Eds.), Proceedings of the Second European Symposium on Biochemical Engineering Science, Univ. of Porto, Porto, Portugal, September 16–19, 93–98, 1998.

[42]. J.C. Maxwell, A treatise on electricity and magnetism, Clarendon Press, Oxford, UK (1873).

[43]. M. Barrande, R. Bouchet, and R. Denoyel, Tortuosity of Porous Particles, Analytical Chemistry, 79, 9115-9121, 2007.

[44]. Q Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Physics of Fluids, 9, 1591-1598, 1997.