Document Type : Full Lenght Research Article
Authors
Semnan University
Abstract
Keywords
Main Subjects
1.
Investigation of porescale random porous media using lattice boltzmann method
Alireza Azhdari Heravi^{*}, Farhad Talebi, Mohammad Sadegh Valipour
Faculty of Mechanical Engineering. Semnan University, Semnan, Iran

History: Received 03 February 2014 Received in revised form 25 September 2014 Accepted 30 December 2014
Keywords: Lattice Boltzmann porescale simulation Creeping flow Random porous media 
A B S T R A C T The permeability and tortuosity of porescale two and threedimensional 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 unitcell as twodimensional porous media, and random granular media in a cubic unitcell as threedimensional 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 twodimensional samples with different porosities and obstacle aspect ratios. Also 29 random spherepackings including three different cases with three different sphere diameters were investigated as threedimensional cases. Employing nonlinear regression based on the “leastsquares” method, two permeability correlations were proposed with minimum curvefitting errors. Besides, the effect of porosity on required timesteps to reach the converged solutions was investigated. It is concluded that an increase in the required timesteps to convergence is seen with reaching both high and low ends of porosity. © 2015 Published by Semnan University Press. All rights reserved. 


Journal of Heat and Mass Transfer Research
Journal homepage: http://jhmtr.journals.semnan.ac.ir 


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 semiempirical formula for the pressure drop which incorporates the tortuosity, the curvature ratio and the variation of the pore crosssectional 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 volumeaveraged 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 latticeBoltzmann 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 threedimensional 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 microscale, 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 porousmedium properties.
In this study, porescale transport properties in the porous medium of various structures are investigated for ﬂows through 2D and 3D porous structures, including 2D unitcell 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 twodimensional 9velocity (D2Q9) and the threedimensional 15velocity (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 postcollision 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 timestep 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 Bounceback BCs :
The halfway bounceback boundary condition is applied on a solid surface [23]. Figure 2 illustrates the halfway bounceback 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 adjacenttoboundary 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 adjacenttoboundary fluid nodes. In this example the bounceback boundary condition is employed as : f_{5=} f_{7, }f_{2=} f_{4}, and f_{6=} f_{8} to make unknown distribution functions for adjacenttoboundary 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
Fig. 2 Halfway bounceback boundary condition [23]. 
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 nonequilibrium 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
Twodimensional domains with a mesh size of 100×100 lu^{2} and threedimensional 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 nondimensionalized by dividing it by the square of the characteristic length scale. For a porescale 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 Twodimensional 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 nondimensionalized 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 opensource simulator.
4.1.1 Permeability
Results by LBM are evaluated against the results of CarmanKozney [27] and modified CarmanKozney 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 deadend 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 CarmanKozney 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) 

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 dashedline 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 deadend 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 Twodimensional dimensionless permeability versus porosity. 

Fig. 9 Effect of obstacle aspect ratios on twodimensional 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 deadend 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 Timesteps
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 timesteps is less than Criteria Number, the simulation will be terminated.
Fig. 10 Twodimensional 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 timesteps 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 timesteps 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 timesteps is increased.
iii. On the other part, figures clearly show that a rise of the number of timesteps is also occurred with a reduction in the porosity. It is because of the fact that at low porosities the volume of deadend pores is increased, which increases the total fluid travel distance and consequently more processing time.
Hence, the timestep is increased along two end limits of the porosity band, and a sharper increase of timestep at high porosities is clearly observed in all three graphs in figure 11.
4.2 Threedimensional Simulations
4.2.1 Permeability
Dimensionless values of extracted permeability in some samples of 100×100×100 lu^{3}spherepacking, versus porosity are shown in figure 12. Results are validated using three prevalent correlations:
i. KozneyCarman [27],
ii. RumpfGupte [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 leastsquares method by fitting the parameters with minimum fitting errors:
(24) 
Table 1. Convergence Number for diverse porosity ranges 

Convergence Number 
Porosity range 
Fig. 11 Timesteps 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 dashline 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 nonporous, nonoverlapping particles [32].
There is another experimental powerlaw correlation which is prevalent for granular media, again with nonoverlapping spheres [3336] with a coefficient of , being dependent on type of spherepacking. According to reports of Refs. [3641], the value of for nonoverlapping 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 nonoverlapping assumption in the experiments. Hence, the present study data, based on overlapping spheres, are fitted to experimental and analytical nonoverlapping formulas, with fitted parameters listed in table 2.
Fig. 12 Threedimensional dimensionless permeability versus porosity. 
5. Conclusions
The flow behavior, permeability, and tortuosity in two and threedimensional random porous media were investigated using Lattice Boltzmann Method with D2Q9 and D3Q15 lattice arrangements. Unitcells with randomly distributed rectangular obstacles, and random spherepacking in cubic unitcells were employed as two and threedimensional 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 twodimensional 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 threedimensional 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 threedimensional 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 twodimensional dimensionless permeability. 
Fig. 14 Threedimensional tortuosity versus porosity. 
Table 2. Fitting parameters of tortuosity. 

Present study, Free to overlap 
Literature, Nonoverlapping 

0.25 
0.4 
Experimental, [32] 
1.4 0.4 
1.5 0.5 
Analytical, [43] 
0.22 
0.4 
Experimental, [3641] 
porosity, but also obstacle aspect ratios. Finally, it is worth mentioning the effect of porosity on the required timesteps 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 opensource 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 LB3DPrime 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;sP2\n100 100\n255\n; s2551g;
s g; s\ng';
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 68.
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, 8994, 1952.
[3]. Macdonald I F, ELSayed M S, Mow K and Dullien F A L 1979 Flow through porous media—the Ergun equation revisited Industrial & Engineering Chemistry Fundamentals, 18 199208, 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), 343360, 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, 3846, 1989.
[8]. S. Kim, W.B. Russel, Modeling of porous media by renormalization of the Stokes equations, Journal of Fluids Mechanics, 154, 269286, 1985.
[9]. A. Nakayama, F. Kuwahara, Y Kawamura, and H. Koyama, Threedimensional numerical simulation of ﬂow through a microscopic porous structure Proc. ASME/JSME Thermal Engineering Conference, 3, 313, (1995).
[10]. D.H. Rothman, Cellularautomaton ﬂ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, 7284, 1991.
[12]. S. Succi, E Foti, and F. Higuera, Threedimensional ﬂows in complex geometries with the lattice Boltzmann method Europhysics Letter, 10, 433438, 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, 20852088, 1990.
[14]. R.S. Maier, D.M. Kroll, Y.E. Kutovsky, H.T. Davis, and R.S. Bernard, simulation of flowthrough bead packs using the lattice boltzmann method, Physics of Fluids, 10 (1), 6074, 1998.
[15]. T. Inamuro, M. Yoshino, and F. Ogino, Lattice Boltzmann simulation of ﬂows in a threedimensional porous structure, International Journal for Numerical Methods in Fluids, 29, 737748, 1999.
[16]. C. Manwart, U. Aaltosalmi, A. Koponen, R. Hilfer, and J. Timonen, LatticeBoltzmann and ﬁnitedifference simulations for the permeability for threedimensional porous media, Physical Review E, 66, 016702, 2002.
[17]. M.A. van der Hoef, R. Beetstra, and J.A.M. Kuipers, LatticeBoltzmann simulations of lowReynoldsnumber flow past mono and bidisperse arrays of spheres: results for the permeability and drag force, Journal of Fluids Mechanic, 528, 233254, 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, 19591978, 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 latticeBoltzmann method, Advanced Water Resources, 29, 13281340, 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 DarcyForchheimer drag for microporous 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 onecomponent systems, Physical Review, 94, 511525, 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.
[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), 406410, 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) , 367375, 1975.
[32]. E. Mauret, M. Renaud, Transport phenomena in multiparticle systems—I. Limits of applicability of capillary model in high voidage bedsapplication to fixed beds of fibers and fluidized beds of spheres, Chemical Engineering Science, 52, 18071817, 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 rocksa 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 Europa97, 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, 91159121, 2007.
[44]. Q Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Physics of Fluids, 9, 15911598, 1997.
[a] Lattice Bhatnagar–Gross–Krook
[b] Or propagation
[c] Hence the pressure by eq. 9
[d] in 2D case, and and in 3D case
[e] Section 2.1
[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.
[i] Which is the unit cell
[j] Created by C++ in binary geometry generation step