Abstract
Discrete element method (DEM) is used to produce dense and fixed porous media with rigid mono spheres. Lattice Boltzmann method (LBM) is adopted to simulate the fluid flow in interval of dense spheres. To simulating the same physical problem, the permeability is obtained with different lattice number. We verify that the permeability is irrelevant to the body force and the media length along flow direction. The relationships between permeability, tortuosity and porosity, and sphere radius are researched, and the results are compared with those reported by other authors. The obtained results indicate that LBM is suited to fluid flow simulation of porous media due to its inherent theoretical advantages. The radius of sphere should have ten lattices at least and the media length along flow direction should be more than twenty radii. The force has no effect on the coefficient of permeability with the limitation of slow fluid flow. For mono spheres porous media sample, the relationship of permeability and porosity agrees well with the K-C equation, and the tortuosity decreases linearly with increasing porosity.
1. Introduction
Fluid flow through porous media has been widespread concern. It plays an important role in the seepage, drainage, and oil engineering. Engineers often need to predict the fluid flow rate of water or oil in porous media. A simple empirical method called Darcy’s law is often used [1]; the formula reads where is outflow speed in the direction , represents volume average, is the viscosity coefficient of the fluid, is the pressure gradient in direction , and is the permeability matrix. Permeability matrix is assumed to be diagonal and positive definite. The elements of permeability matrix can be measured by experimental methods or estimated by computational methods.
Many factors affect permeability, such as porous media structure, fluid physical properties, and degree of fluid saturation. This paper studies the impact of porous media structure only. The traditional predictive research is usually concerned with the relationship between permeability and porosity. In addition to that, specific surface of the particle, shape of the particle, and tortuosity also have a significant impact on permeability. But effects of these factors are hard or even impossible to measure accurately; these problems need pore-scale simulation of fluid flow in porous media.
Because of the complexity of fluid flow in porous media, except for some very simple question, it is difficult to get the analytical solution. With the development of the porous media geometry visualization, computing power, and computational methods, it is possible to simulate directly fluid flow in pore-scale porous media.
The LBM method has been applied to this field since it was introduced in 1989. Ghassemi and Pak [2] studied the effects of permeability and tortuosity on flow through particulate media by LBM and DEM; the relationships between permeability and tortuosity with other parameters such as particles diameter, grain specific surface, and porosity were identified. The results show that the relationship between permeability and porosity of the porous media agreed the K-C formula well, but there is a bit deviation in the high porosity region; tortuosity almost linearly decreases with the porosity, and the particle diameter does not affect tortuosity any longer after it becomes large enough.
Sun et al. [3] introduced a multiscale method. Microstructural attributes such as occluded/connected porosity and geometrical tortuosity are extracted using new computational techniques to form digital images of porous media, and a LBM/FEM scheme is used to obtain homogenized effective permeability at specimen scale. The results show that the permeability depends not only on the porosity, but also on the connectivity and the tortuosity.
Experimental measurements and LBM simulation on samples of packed particle beds were researched by Videla et al. [4]. Due to the resolution limit, when the particle diameter is greater than 500 microns, the experimental results and numerical results are consistent. The impact of the relaxation time on permeability was also studied. Vidal et al. [5] verified that the LBM numerical permeability of the porous media of spherical particles agrees with experimental data very well. When the size distribution of spherical particles disperses to a certain extent, the results start to deviate from K-C formula, and a modified formula was given according to the results. Jeong [6] studied the micro flows through granular porous media at various Knudsen numbers, the correlation between the permeability, the porosity, and the Knudsen number is derived, a new fitting formula is given, and the effect of rarefaction on the permeability is also discussed. van der Hoef et al. [7] simulated fluid flow past mono- and bidisperse random arrays of spheres with LBM, the calculated permeability is consistent with experimental results, and a new drag force formula is given. The effect of the arrangement of the particles and particle shape on permeability anisotropy was researched by Stewart et al. [8].
In this work, the permeability and tortuosity of fluid flow in three-dimensional dense porous media have been simulated and analyzed. The effects of porosity and the size of particle have been studied. The DEM is described in Section 2. Numerical model of LBM is present in Section 3. Calculation results are given in Section 4 and the Conclusion in Section 5.
2. Discrete Element Method and Porous Media Generation
2.1. Discrete Element Method (DEM)
DEM is an explicit numerical method for modeling motion of distinct particles. Since this method was introduced by Cundall and Strack [9], it has become a powerful tool and has been used to research a wide range of problems in science and engineering.
In DEM, the granular porous media are regarded as an assembly of distinct rigid particles. When overlap occurs between particles (or particle and wall), the contact force separates them from each other. Using the force-displacement law, the force system of each particle can be summed. The displacement and velocity of each particle can be calculated by Newton’s second law of motion.
The force-displacement law gives the relationship between the contact force and the overlap between particles (or particle and wall) by the following equations: where , , and are increments of force, stiffness, and overlap between particles, respectively. The superscripts and represent normal and tangential direction of overlap plane.
The contact forces consist of normal force and tangential force. The tangential force is related to the normal force by the Coulomb frictional law. If , relative translation motion will not occur. Here is the sliding friction coefficient between the contact particles, and , are the normal and tangential forces on particle in position .
Newton’s second law of motion describes the relationship between the acting forces and movement of particle. For sphere , the law reads where , are the mass and inertia of sphere particle ; and are linear and angular accelerations in direction ; , are the global damp force and moment; is the vector from the center of sphere to the contact point.
2.2. Porous Media Generation Using DEM
The first step of modeling is producing dense assembly of granular particles. We retrieve it with opening LAMMPS code. The details are shown in Figures 1 and 2. A domain whose size is by lattice unit is created first. Periodic boundary condition is used in the two horizontal directions. In vertical direction, rigid and impassable planes are set at the top and bottom boundaries.
The rigid particles are created in a small region near the top boundary first, in order to reduce the probability of overlap; the ratio of particles volume to domain is set to ten percent. The position of each particle will be determined when it does not overlap with any other particles. Because the percentage of particles is low, this process is easy to implement.
After the loose assembly is created in the specific region, it starts to fall down under gravity (Figure 1). When the particles leave the initial region totally, new particles will be created. This process is repeated until enough particles appear in computational domain. The particles interact with other particles, and the dense assembly of particles forms under gravity at last.
To get loose porous media with a desired porosity, the radii of particles are shortened with same scale, but the amount and positions of particles are unchanged.
3. Lattice Boltzmann Method and Simulation of Fluid Flow
3.1. LBM
The simplest LBM model is known as the BGK-Boltzmann model. It was first introduced by Bhatnagar et al. [10]. The evolutional equation reads where , , , and are probability distribution function of particle in direction , space position, time, and microscopic velocity, respectively. moves from to in a time step. The right term indicates the colloids function, that is, relaxation function. Probability distribution function of particle approaches (the local equilibrium distribution function). is the function of local density and velocity. The relaxation factor indicates relaxation speed rate and is related to fluid viscosity coefficient by . Every evolution procedure consists of two stages, a propagation stage and a collision stage. Under limitation of low fluid flow, the steady solution of (4) is the solution of Navier-Stokes equations.
The final fluid flow state is evidently influenced by boundary conditions. The interface of the particles and fluid is assumed as a no-slip boundary. The simplest method which is called half-way bounce-back method is adopted here. The pressure drop is implemented by placing a body force on fluid. Periodic boundary condition is used in all directions. The D3Q19 model is used to calculate the fluid flow among the interval of porous media (Figure 3). This model includes a stationary part and eighteen moving parts. They are ; here and .
Local equilibrium distribution function is given by where , , , and are density, sound velocity, velocity vector, and constant coefficient, respectively. The sound velocity equals and the constant coefficients are
New density and velocity of fluid can be extracted from probability distribution function of particle by equations
3.2. Simulation of Fluid Flow Using LBM
As input condition, assembly of dense particles created by DEM is transmitted into LBM code. Lattice is signed by TRUE if it does not locate inside any particle, FAULSE if it does. All particles are fixed during simulation procedure, and they interact with fluid by bounce-back boundary condition [11].
Lattice should be fine enough to catch the details of the particles. We place ten lattices along the radius of sphere in our simulation. The computer load is quite big, and an elementary lattice structure is .
In flow direction, additional layers are added to the ends of domain. The path of fluid flow is shown in Figure 4. Periodic boundary condition is used in the other two directions.
A uniform body force drives fluid to flow; the GZS model is adopted in our simulation [11]. All the simulations start from a zero velocity filed; convergence to the steady solution is examined by the following criterion: The summation covers all the fluid lattice sites, and superscript denotes the time steps.
4. Results and Discussions
The permeability of the porous media can be computed by formula (1) from steady flow of fluid. The tortuosity is defined by , where and are the velocity vector and the velocity component in flow direction at position , respectively.
Permeability and Lattice Number of the Radii. The space is dispersed into regular lattice point in the LBM. The porous media particles are described by the lattice point. In order to accurately capture the shape of the particles, enough lattices must be placed along the radius of the particle. The more the lattices are, the more accurate the particles are depicted, but it also brings the problem of excessive calculation. Therefore, we need to find a balance point between accuracy and computation burden.
In this paper, the lattice unit is applied always. To simulating the same physical problem, the key criteria numbers of simulations and real physical problem should be matched. These key criteria numbers are reference length , reference density , and reference velocity . They are read here , , , and are length, density, time, and viscosity in lattice unit respectively, and , , , , and are correspond parameters in physical unit respectively.
Therefore, the viscosity and body force must be adjusted simultaneously. The permeability is obtained with different lattice number. Figure 5 shows the relationship between permeability and the number of lattices among the radius; the permeability enlarges with increase of lattice number. For controlling the computational time (considering machine memory and convergence time), ten lattices are placed along the radial direction (red data points). According to Videla et al.’s [4] work, the difference between experimental measurements and numerical simulations is less than 30% when the ratio of diameter to resolution is bigger than 25. If the ratio is less than 20 (about 17.5 and 12.5), the difference is very big. Our numerical results are consistent with Videla’s conclusion.
Permeability and Media Length. Along the flow direction, the sample of porous media should guarantee enough length. The permeability of the media should be independent of media length in the flow direction. Fluid flow in porous media for different lengths of sample is simulated (see Figure 2); Figure 6 shows the variations of permeability versus the media length in the flow direction. From the figure, when the ratio of channel length to radius is more than 20, variation of permeability is very small. This media length is used in all our simulations (red data points).
Permeability and Body Force. Similarly, the permeability of porous media should not depend on the body force (pressure gradient). Figure 7 depicts the permeability results of the same sample under different body force. The permeability remains unchanged when the body force is small; the permeability is significantly increasing when the body force becomes larger. It should be noted that maximum speed is 0.0124 when the body force is . In the following calculation, the body force is taken as (the red data point in Figure 7).
Permeability and Porosity. About relationship between porosity (void ratio) and permeability, researchers proposed some different equations. Carman introduced and extended the K-C equation [12]; the equation reads is a constant parameter which is related to the specific surface area, tortuosity, and fluid physics, and is void ratio. Ghassemi’s simulation gave [2]. In order to study the influence of porosity on permeability, a series of porous media samples with uniform sphere particles is simulated by LBM. All the samples are made of sphere particles with the same diameter, but they are different in porosity. In Figure 8, the relationship between regularization of the permeability and void radio is shown. The curve fits well with the K-C equation. Using least squares method, we have . When the void ratio is less than 1, the three results are similar; for large porosity, our numerical results are also consistent with K-C equation, but are different from Ghassemi’s results. Ghassemi’s simulation is based on two-dimensional porous media, and his model cannot capture all the physical characteristics of real porous media; the author elaborated this fact in his paper.
Permeability and Particle Size. Figure 9 depicts the relationship between permeability and particle size; the porosity remains constant () in all samples. According to the K-C equation, the permeability is inversely proportional to the square of the specific surface area. For the homogeneous sphere porous media, the permeability should be proportional to the square of the diameter. From our results, the relationship between the permeability and the radius does not seem to comply with this rule. In fact, it is because the range of particle radius is narrow and the permeability is basically proportional to the radius.
Tortuosity and Porosity. For the relationship between tortuosity and porosity, different researchers have given different results also. Koponen and Ghassemi’s results show that the tortuosity decreases linearly when porosity increases [13]. Yet, Boudreau gave the nonlinear relationship [14]; Dias believes that exponential relationship is better [15]. For porous media of uniform particles, Figure 10 depicts the relationship between tortuosity and porosity. The tortuosity essentially decreases linearly with increasing porosity. But our result is much lower than the absolute value of Koponen and Ghassemi’s results. The reason is that the dead end will not form in the three-dimensional porous media samples. For example, when the porosity is equal to 0.5, Koponen and of Ghassemi’s results are 1.4 and about 1.33, respectively; our result is 1.14. For the slope of the curve, the results of our calculations (about −0.42) and Ghassemi (−0.45) are relatively close and far from Koponen’s result (−0.8). The porous media in Koponen’s paper are composed by small rectangles. The channel among porous media has more twists and turns, and the permeability is more sensitive to the porosity. Therefore, when the porosity increases, the tortuosity decreases faster.
5. Conclusion
DEM and LBM are used to simulate the fluid flow in the dense porous media; the ability of LBM to resolve this problem is confirmed. The relationship between permeability, tortuosity and porosity, and particle size is researched. In our calculation, all the sample particles are uniform size, rigid, and fixed; high porosity samples are created through rescaling of dense samples. We obtain the following conclusions.
LBM is suitable for the simulation of flow in porous media because of its inherent advantages. The results show that the radius of the sphere should not be less than ten lattices; the ratio of channel length to radius should be larger than 20. Under the limitation of low speed, body force has little influence on the calculation results. If relaxation time is equal to one, the convergence can be achieved quickly, or the convergence is extremely slow. For porous media with uniform spheres, the relationship of permeability and porosity basically meets the K-C formula, and the tortuosity decreases linearly with increasing porosity.
Simulation of samples with a wide range of particle size will be run in further work, and the permeability calculation of nonuniform and nonspherical particle samples will be studied also.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The project was supported by the National Natural Science Foundation (11102034) and the Important National Science & Technology Specific Project (2011ZX02403-004).