Azhdari Heravi, A., Talebi, F., Valipour, M. (2015). Investigation of pore-scale random porous media using lattice boltzmann method. Journal of Heat and Mass Transfer Research(JHMTR), 2(1), 1-12. doi: 10.22075/jhmtr.2015.326

Alireza Azhdari Heravi; Farhad Talebi; Mohmmad Sadegh Valipour. "Investigation of pore-scale random porous media using lattice boltzmann method". Journal of Heat and Mass Transfer Research(JHMTR), 2, 1, 2015, 1-12. doi: 10.22075/jhmtr.2015.326

Azhdari Heravi, A., Talebi, F., Valipour, M. (2015). 'Investigation of pore-scale random porous media using lattice boltzmann method', Journal of Heat and Mass Transfer Research(JHMTR), 2(1), pp. 1-12. doi: 10.22075/jhmtr.2015.326

Azhdari Heravi, A., Talebi, F., Valipour, M. Investigation of pore-scale random porous media using lattice boltzmann method. Journal of Heat and Mass Transfer Research(JHMTR), 2015; 2(1): 1-12. doi: 10.22075/jhmtr.2015.326

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

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

Alireza Azhdari Heravi^{*}, Farhad Talebi, Mohammad Sadegh Valipour

Faculty of Mechanical Engineering. Semnan University, Semnan, Iran

PAPER INFO

History:

Received 03 February 2014

Received in revised form 25 September 2014

Accepted 30 December 2014

Keywords:

Lattice Boltzmann pore-scale simulation

Creeping flow

Random porous media

A B S T R A C T

The permeability and tortuosity of pore-scale two and three-dimensional random porous media were calculated using the Lattice Boltzmann method (LBM). Effects of geometrical parameters of medium on permeability and tortuosity were investigated as well. Two major models of random porous media were reconstructed by computerized tomography method: Randomly distributed rectangular obstacles in a unit-cell as two-dimensional porous media, and random granular media in a cubic unit-cell as three-dimensional porous media. Results were validated using available theoretical, experimental, and numerical results from the literature. It is observed that permeability is a weak function of porosity in low porosity regions, but a strong function of porosity at high porosities. It also depends on the aspect ratio and hydraulic diameter of obstacles.

Permeability results were obtained regarding to 73 random two-dimensional samples with different porosities and obstacle aspect ratios. Also 29 random sphere-packings including three different cases with three different sphere diameters were investigated as three-dimensional cases. Employing nonlinear regression based on the “least-squares” method, two permeability correlations were proposed with minimum curve-fitting errors. Besides, the effect of porosity on required time-steps to reach the converged solutions was investigated. It is concluded that an increase in the required time-steps to convergence is seen with reaching both high and low ends of porosity.

Journal of Heat and Mass Transfer Research 2 (2015) 1-12

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 modiﬁcations 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 ﬂows, Larson and Higdon [7] performed calculations for the Stokes ﬂow 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 ﬂows 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 ﬂow passages. In recent years, the Lattice Boltzmann Method (LBM), as a mesoscale approach, has emerged as an alternative tool to investigate ﬂow in complex geometries. In general, the Lattice Boltzmann Method (introduced in the next section) is easier to be implemented than conventional computational ﬂuid dynamic techniques, is highly compatible with parallelization, and can deal with arbitrary complex ﬂow geometries without heavy penalties. These attributes, combined with its ﬂexibility, make LBM suitable for the investigation of a wide range of ﬂow problems ranging from flows through complex geometries (e.g. permeating groundwater ﬂow, melt segregation) to ﬂows in which multiphase interactions are of interest (e.g. particle laden ﬂows, 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 ﬂows 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 ﬂows 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 ﬂow regime even though the LBM can be applied to the slip ﬂow 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 propertiesin the porous medium of various structures are investigated for ﬂows 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/lbﬂow/. 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

,

(1)

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

(2)

(3)

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

(4)

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

Fig. 1 D2Q9 (left) and D3Q15 (right) lattice velocity directions

(5)

(6)

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

The collision step is:

(7)

and streaming step is :

(8)

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

(9)

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]:

(10)

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_{2}, f_{5} and f_{6} 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_{5=} f_{7, }f_{2=} f_{4}, and f_{6=} f_{8} 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, :

Fig. 3 A schematic illustration of pressure boundary condition at the inlet boundary (density is known).

Note that velocity tangent to the boundary, u_{y}is 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:

(11)

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

(12)

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:

(13)

where and denote inlet and outlet, respectively.

3.2 Geometry Generation

Two-dimensional domains with a mesh size of 100×100 lu^{2}and three-dimensional domains with a 100×100×100 lu^{2} mesh size were selected with periodic boundary conditions in all three spatial directions. To induce the ﬂow, 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

Fig. 4 A schematic illustration of the periodic boundary condition.

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.

(14)

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

Now the dimensionless permeability is defined as :

(15)

Fig. 5 Binary Geometries; Left: 3D, Right: 2D.

Fig. 6 Part of an ASCII geometry file.

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:

(16)

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.

(17)

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:

(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,

(19)

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.

(20)

Fig. 7 Contours of velocity magnitude in 2D random geometries with AR=0.5(left), 1(middle), and 2(right).

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]:

(21)

Fig. 8 Two-dimensional dimensionless permeability versus porosity.

Fig. 9 Effect of obstacle aspect ratios on two-dimensional dimensionless permeability.

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:

(22)

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

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

(23)

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.

Fig. 10 Two-dimensional tortuosity versus porosity.

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^{3}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:

(24)

Table 1. Convergence Number for diverse porosity ranges

Convergence Number

Porosity range

Fig. 11 Time-steps required for macroscopic local velocities to converge.

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.

Fig. 12 Three-dimensional dimensionless permeability versus porosity.

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

Fig. 13 Effect of obstacles’ diameter on two-dimensional dimensionless permeability.

Fig. 14 Three-dimensional tortuosity versus porosity.

Table 2. Fitting parameters of tortuosity.

Present study, Free to overlap

Literature, Non-overlapping

-0.25

-0.4

Experimental, [32]

1.4

0.4

1.5

0.5

Analytical, [43]

0.22

0.4

Experimental, [36-41]

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

Fig. A1 Creating 2D binary geometries with Microsoft Paint

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|P2n100 100n255n||; 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.

Fig. A2 (repeated). Part of a resulted ASCII file

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 ﬂow 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 ﬂow of ﬂuids 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 ﬂow 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 ﬂow through a microscopic porous structure Proc. ASME/JSME Thermal Engineering Conference, 3, 313, (1995).

[10]. D.H. Rothman, Cellular-automaton ﬂuids: a model for ﬂow 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 ﬂow through porous media, Physica D, 47, 72-84, 1991.

[12]. S. Succi, E Foti, and F. Higuera, Three-dimensional ﬂows 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 ﬂows 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 ﬁnite-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, LBﬂow: An extensible lattice Boltzmann framework for the simulation of geophysical ﬂows. Part I: theory and implementation”. Computers & Geosciences, 36, 115–122, 2010.

[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.

[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.

[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.

[f] Thus we can consider forchheimer drag, F.Re_{D}, negligible.

[g] A comprehensive explanation for this correlation can be found in [43].

[h] The software measures the total area occupied by black and white portions, separately. Accordingly, the ratio of black area to the total area, porosity, is calculated easily by the software itself.

[j] Created by C++ in binary geometry generation step

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 ﬂow 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 ﬂow of ﬂuids 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 ﬂow 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 ﬂow through a microscopic porous structure Proc. ASME/JSME Thermal Engineering Conference, 3, 313, (1995).

[10]. D.H. Rothman, Cellular-automaton ﬂuids: a model for ﬂow 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 ﬂow through porous media, Physica D, 47, 72-84, 1991.

[12]. S. Succi, E Foti, and F. Higuera, Three-dimensional ﬂows 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 ﬂows 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 ﬁnite-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, LBﬂow: An extensible lattice Boltzmann framework for the simulation of geophysical ﬂows. Part I: theory and implementation”. Computers & Geosciences, 36, 115–122, 2010.

[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.

[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.

[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.