Abstract
The hydraulic conductivity () represents an important hydrophysical parameter in a porous media. direct measurements, usually demand a lot of work, are expensive and time consuming. Factors such as the media spatial variability, sample size, measurement method, and changes in the sample throughout the experiment directly affect evaluations. One alternative to measurement is computer simulation using the Lattice Boltzmann method (LBM), which can help to minimize problems such as changes in the sample structure during experimental measurements. This work presents experimental and theoretical results (simulated) for three regular finite arrangements of spheres. Experimental measurements were carried out aiming at corroborating the LBM potential to predict once the smallest relative deviation between experimental and simulated results was 1.4%.
1. Introduction
Hydraulic conductivity () is an important parameter in processes of fluid flow in porous media. indicates how easily certain fluid is transported through porous media, and depends on the media properties as well as percolating fluid characteristics. Pore size distribution, type of pores, tortuosity, and connectivity are some of the factors related to the porous media. Regarding the percolating fluid, its viscosity () is the main factor related to measurements [1]. For example, increase in water temperature reduces its viscosity and potentially increases .
determinations are usually characterized by great variability due to factors such as media spatial variability, sample size, measurement method, changes in the sample throughout the experiment, and others [2]. Therefore, representative determination of requires several measurements and samples. Direct measurements are usually expensive and demand hard and thorough technical work [3].
Theoretical models and numeric simulations which enable measurement from information about the porous media structure might be an interesting alternative to predict this physical parameter [4–7].
One theoretical tool that can be successfully used to predict is the Lattice Boltzmann method (LBM) [8, 9]. The LBM is based on evolution of a relaxation equation for fluid particles distribution function, which is related to density and fluid macroscopic momentum. In the LBM equation, there is input data, which is the relaxation time that is related to the number of time steps so that the thermodynamic equilibrium is reached defining fluid viscosity [10]. The LBM can reproduce the macroscopic behavior of a fluid according to the Navier-Stokes (NS) equation. Relative easiness of computation implementation and numeric stability, in a great variety of flow conditions, makes the LBM ideal for treatment of fluid flow in porous media [11].
The development of LB models had important advances in recent decades; for example, today it is possible to simulate compressible, heat, and multiphase flows [10, 12]. However reliability of these applications has occurred through validations without regard to analytical results, which are usually associated with simplified cases of reality. There are few studies in which the LBM is validated with experimental results and the structure of porous media to be simulated is usually obtained indirectly by techniques (computed tomography or image reconstruction) that have their intrinsic sources of error [9, 11, 13].
LBM has been used for problems in porous media under several aspects [9, 11, 14, 15]. However, such studies have posed extra difficulties to the LBM such as image capture and reconstruction of the media as a representative porous media and results presented take into consideration possible deviations related to these difficulties.
In order to create higher possibility of comparison between results of LBM simulations and experimental results, this work suggests simulating the fluid flow through a layer of spheres. This is due to the fact that spheres have the simplest form to be digitally built and their symmetry enables the control of their superficial irregularity on the results, since by increasing their diameter such irregularities are less perceived by the flow.
We propose an experiment where the construction of the porous media is greatly facilitated—a finite array of spheres. The unique source of error is the roughness (discretized surface) of the sphere, which can be controlled with its diameter increase. So, we present in this work experimental and simulated results of measurements in three porous media constituted of regular arrangements of spheres.
The main objective is to show that the LBM method can be used to evaluate in the arrangements analyzed. This objective was achieved through comparisons between experimental and simulated . The success of results proposed in this work points to the future use of this method in representative measurements of more complex porous media such as soil samples.
2. Materials and Methods
2.1. Experimental Methods
An experimental apparatus tested by Camargo et al. [8] was used for K experimental measurements (Figure 1(a)). The steps below were followed to determine : (a) porous media saturation (acrylic box with a certain sphere arrangement) with glycerin, C3H5(OH)3, (manufactured by Biotec, 99.5% purity); (b) length measurement (hydraulic load) and (porous media height); (c) percolated glycerin mass measurement to obtain its volume, where the glycerin density is known; (d) measurement of the necessary time interval for the glycerin to percolate; (e) use of Darcy’s Law to calculate conductivity using = volume of percolated glycerin (cm3), = cross sectional area of the box containing spheres (cm2), and = time interval for a given volume of glycerin to percolate (s) [16]; (f) glycerin viscosity measurement using its flow through a D = 0.27 cm diameter acrylic cylinder and the analytical expression for the cylinder conductivity : = gravity acceleration (cm s−2); (g) was measured using the steps (a) to (e).
(a)
(b)
(c)
In both measurement cases a sufficiently big was guaranteed so that the glycerin would not form drops when leaving the spheres for the measurement, or the cylinder for the viscosity measurement, which would delay the measurement time used in Darcy’s equation.
Glycerin was used because it presents high ν leading to a Reynolds number (Re) smaller than 1, where Darcy’s Law is valid [17]. Once the glycerin is highly susceptible to temperature variations, its measurement was carried out for each layer of spheres added during the experimental arrangements (Figure 1(b)).
2.2. The Lattice Boltzmann Method
In order to simulate through the LBM a 3D media was built (Figures 2(a)–2(c)) similar to the real porous media being represented in a binary language 0 (solid) and 1 (porous) distributed along the vertices of a regular lattice. Once the porous media was built, the computer program that simulates (2) was used. The program returns the media’s intrinsic permeability () as well as the pressure and flow velocity fields (Figure 1(c)).
(a)
(b)
(c)
(d)
The lattice used in the simulations was the cubic D3Q19. The 18 direction vectors of this lattice (Figure 2(d)) connect the sites one to another and also represent the possible velocity vectors, and there is still the null velocity (19 velocities).
Being a site in the lattice located by the vector and having close neighbors, the evolution equation for the particle of fluid distribution function is given by the Lattice Boltzmann Equation: = vector coordinates of the site in the lattice (lattice units), = time step variable . The duration of a time step is taken to be unity that represents the time interval for the particle of fluid travel between the closest neighbors, = number of fluid particles (direction ) located at site at time , = collision operator that represents the collision of fluid particles with others at time (see (4)), = direction of one of the closest neighbors , and = velocity vector in direction . This vector coincides with the lattice vectors, because in a unit time step, a particle travels from one site to adjacent one. The term represents the resting particles.
Variables and are given in the called lattice units and scale factors are necessaries for these variables assume length () and time (δ) dimensions, that is, it can be assumed that, for example, 5 units of the lattice is equivalent to 1 mm (h = 1 mm/5) or 5 time steps equivalent to 1 second (δ = 1 s/5). Thus, the scale factor for the velocity vector becomes h/δ.
In this work we assume that variables without units will be represented as lattice units, that is, length and time variables have as unit the lattice spacing and the time step, respectively. So, velocity, viscosity, pressure, and other properties will be represented by lattice units.
The mesoscopic dynamics occurs in two steps: (1) propagation step represented by (2); (2) collision step represented by (3), which simulates the molecular collisions needed so that thermodynamic equilibrium occurs. This step is given by the action of collision operator on the : = “collided” distribution function that will present a new value (number of fluid particles) at site , in direction, and time .
A simple and sufficient form of collision operator which recovers the Navier-Stokes macroscopic equation is known as BGK (variables and were omitted from here to reduce notations) operator [19]: = relaxation time, which is a function of fluid viscosity, = equilibrium distribution (see (8)).
So, if , and the amount will be added to making tend to .
The macroscopic particle density and the macroscopic momentum in a site are given by: = density at site at time .
Taking that into consideration, the collision operator conserves mass and momentum, where .
Particle distribution for the is usually obtained through the expansion, in power series at the macroscopic velocity (), being sufficient so that Navier-Stokes equation is recovered. For the low Mach number the pressure () is given by: where is the Euclidian dimension of space in which the lattice is immerse and . With this, the balance distribution form for moving particles is given by:
In the main directions , , and , the balance distribution must be doubled () so that viscosity is isotropic [20].
Resting particles have the following balance distribution:
For a macroscopic analysis of the dynamics proposed by (2), time and space scales are usually used and the Knudsen variable is defined, where and are, respectively, the macroscopic characteristic length and time. With this, the Chapman-Enskog method [20] can be used, considering the equilibrium distribution disturbance, to show that (2) becomes the Navier-Stokes, given by (10), disregarding the O() contributions, in which and are indexes which represent spatial coordinates x, y, or z; for these indexes Einstein’s notation is seen (sum over repeated indexes). Equation (10) is the component of the Navier-Stokes equation with kinematic viscosity given by:
So when the hydrodynamic limit is imposed , the macroscopic dynamic given by (10) is reproduced by (2).
The boundary conditions used in the simulation are periodic, that is, the fluid which leaves one end of the cavity is injected in the other end. The interaction between fluid and solid occurs so that there is no sliding, in this case, the “bounce back” condition was adopted, where the fluid which collides with the walls has its velocity inverted. The program returns permeability (), which is calculated (12) with the Santos method [11], where in a stationary flow, the strength applied to the fluid is equal to the loss of momentum on the walls: = porosity, = fluid average momentum in the porous media, = fluid average momentum lost in the collisions with the porous media walls in a propagation step.
In order to obtain from , the relation below is used: = permeability, = kinematic viscosity (see (11)).
2.3. Computational Simulations
Different diameters (D) of spheres were investigated in the simulations (Figure 3), in order to better compare theoretical results with the ones obtained experimentally (Figure 4).
(a)
(b)
(c)
(d)
(a)
(b)
(c)
The contour conditions used in simulations are periodic and make the fluid that leaves one end of the simulation dominium to enter the other end, as in an infinite array of spheres. Therefore, in the QQ arrangement simulation (Figure 2(a)) the k permeability of a single layer of spheres is the same as the permeability of any other number of layers. Despite that, the conductivity might be different, once viscosity is altered by the room temperature in the experimental measurement of each layer.
3. Comparison between Experimental and Simulated Results
Larger D (Figure 4) provides better approximation between simulated and experimental results. This is due to the reduction in the sphere surface discretization. In the QQ arrangement (Figure 4(b)), porosity is kept constant by the arrangement symmetry, in the QU (Figure 4(c)) the same does not occur and it is necessary to simulate the flow on all sphere layers, which was carried out for D = 34 (scale factor h = 0.3175 mm/34). Concerning the OO arrangement (Figure 4(a)), the spheres have twice the diameter of the two previous cases. The discrepancy from the first layer to the others remains and is a case to be further investigated.
Experimental hydraulic conductivities were plotted against the computed based on the LBM simulations (Figure 5). The data obtained clearly shows that based on computed simulations are in very good agreement (high positive correlation coefficients) with the laboratory measurements. Analyzing the correlation graphics (Figures 5(a) and 5(b)) most of the results are within the variability limits of the laboratory measurements as indicated by the experimental error bars.
(a)
(b)
In relation to discrepancies found in conductivity values from the first to the other layers, simulations were carried out with a buffer zone up to 100 sites length at the input and output of the flow to check the influence of the periodic contour condition, which does not exist in this experiment. However, no alteration was detected in the conductivities. The interfacial tension might have some effect in the flow output since in the cases in which the hydraulic load is not big enough there is some glycerin dropping. When the hydraulic load is big enough, a stream forms in the flow output, where there are two fluids, glycerin and air; the same does not happen in simulations where there is a monophasic flow. Although there are collision operator models to simulate biphasic flows [10], it is still complicated to control the huge difference between the viscosities of glycerin and air.
Due to some limitation of the PC RAM memory recognition by the software used (~2 Gbytes), it was only possible to simulate one layer of spheres with a maximum diameter D = 122 (the same first layer as the QQ case). However, the maximum and minimum deviation for D = 34 were 18.6% and 15.5% for QU. It is relevant to point out that the deviations obtained might be minimized with the use of clusters and parallel processing for simulations with larger diameter spheres, mainly in the QU case which represents a more complex media than QQ and OO.
4. Concluding Remarks
In this work, an experimental physical reality was pursued with flow in Reynolds low number and spherical symmetry to facilitate the media digital construction, avoiding image acquisition problems. This provided a close comparison between the experimental measurement and the one simulated via LBM, with maximum and minimum deviation (using the experimental value as a reference) from 18.8% and 1.4% for QQ, with the maximum deviation happening only in the first layer. For the OO case, the deviations were 35.6% and 3.9% with the maximum deviation being observed in the first layer again.
In the soil science area, LBM can be used in association with the X-ray computed tomography (CT) utilized to acquire more real 3D soil pore structures. The selection of adequate image analysis procedures, for example, threshold, will allow to accurately reconstructing the pore system structure used to simulate for heterogeneous and nonsymmetrical media such as soil. It is important to mention that no extra computational difficulty is included in this case.
The use of 3D soil images will make it possible to simulate the 3D fluid flow allowing the evaluation of important hydraulic soil properties such as and . As direct measurements, usually demand a lot of work, are expensive and time consuming, LBM can be an interesting tool to its simulation. With LBM it will also be possible to access the computed flow velocity, which can be utilized for instance to better understand some important phenomena that occur in the soil such as water fingering.
Acknowledgments
The authors would like to thank the Brazilian agencies Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) for the research fellowship to L. F. Pires and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) for the M.S. fellowship granted to M. A. Camargo.