Simulation of Non-Newtonian Fluids Through Porous Media

The past two decades have appeared techniques for simulation of hydrodynamic systems (Mc Namara & Zanneti, 1998; Higuera & Jiménez, 1989; Qian et al., 1992; Doolen, 1990; Benzi et al., 1992), magnetohidrodynamics (Chen, 1991), multiphase and multicomponent fluids (Shan & Chen, 1993), included suspensions (Ladd, 1994) and emulsions (Boghosian et al., 1996), flows chemically reactive (Chen et al., 1995) and multicomponent flows through porous media (Landaeta, 1997), only for mention some. This techniques, don’t based on directly in the discreet state of hydrodynamic equations, don’t ether tally with the microscopic level of molecular dynamic. These techniques have been named techniques of mesoscopic simulations. In general, these discreet techniques involve collisions of “particles” which retain mass, momentum and some cases energy and consistently, give place a hydrodynamic macroscopic behavior. These methods are known like Cellular Automaton and Lattice Boltzmann.


Introduction
The past two decades have appeared techniques for simulation of hydrodynamic systems (Mc Namara & Zanneti, 1998;Higuera & Jiménez, 1989;Qian et al., 1992;Doolen, 1990;Benzi et al., 1992), magnetohidrodynamics (Chen, 1991), multiphase and multicomponent fluids (Shan & Chen, 1993), included suspensions (Ladd, 1994) and emulsions (Boghosian et al., 1996), flows chemically reactive (Chen et al., 1995) and multicomponent flows through porous media (Landaeta, 1997), only for mention some.This techniques, don't based on directly in the discreet state of hydrodynamic equations, don't ether tally with the microscopic level of molecular dynamic.These techniques have been named techniques of mesoscopic simulations.In general, these discreet techniques involve collisions of "particles" which retain mass, momentum and some cases energy and consistently, give place a hydrodynamic macroscopic behavior.These methods are known like Cellular Automaton and Lattice Boltzmann.
The Lattice Boltzmann Method (LBM), including the method Cellular Automaton (AC), present a powerful alternative to standard approaches known like "of up toward down" and "of down toward up".The first approximation study a continuous description of macroscopic phenomenon given for a partial differential equation (an example of this, is the Navier-Stokes equation used for flow of incompressible fluids); some numerical techniques like finite difference and the finite element, they are used for the transformation of continuous description to discreet it permits solve numerically equations in the computer.
The second approximations study a microscopic description of the particles, through molecular dynamic's equations.Here the position and speed of each atom or molecule in the system are calculated with the solution of the Newton's movement equations.
Between the two approximations, are the LBM and the AC, which are considered mesoscopic approaches that was mentioned previously (Wolf-Gladrow, 2000).
With regard to the LBM, this has his origin en the AC, the first effort for simulate systems of flow of fluids was made in 1973 by Hardy, Pomeau and de Pazzis (HPP) (Hardy et al, 1973), whom showed the model into and this is named with their initials "HPP".In this model, all particles have mass unit and speed; the particles are limited to travel in directions With regard to the LBM, it has its basis in concepts of the kinetics of gases theory.Is method has been employed for simulate many physical phenomenon with the use of discrete lattice, the Boltzmann equation and computational tools.For example in dynamic of fluids, the Navier-Stokes equations can't solve directly, but if when recover like consequence.The problem is obtaining the function of particle's distribution, ( )  i fx , t and starting from these obtains observable variables like viscosity, fallen pressure and Reynolds number.In general, the methods of lattice Boltzmann consist in two operations; the first is denoting the advance of particles at the neighbouring site of the lattice.The second operation, is for simulate particle's collisions.An important advantage is that this method is fully parallel and local; too is simple made the programming.In this can incorporate boundary conditions and complexes geometries (Chen et al, 1994).

Basis of kinetics theory and Boltzmann equation
Fort the solution of complex problems of dynamics of fluids, exist traditionally two kinds of points of view: the first is macroscopic, which is considered continuous, with an approach of differential equations in partial derivatives, for example of Navier-Stokes equations used for flow of incompressible fluids and numerical techniques for its solution.The second point of view is microscopic it has its basis in kinetics theory of gases and statistical mechanics.
With regard to the kinetics theory of gases, this mentions that a gas contain many particles in interaction α ( t h e o r d e r i s 1 0 23 ), the physical state is describing by their positions r, v and each particle can be describing through its trajectory in the space of 6-dimensional by complete extension α  r and α  v , it is named phase space (Struchtrup, 2005).
This chapter will shortly explain the basics of the Boltzmann equation, the derivation of the Lattice Boltzmann equation and the connection to the Navier-Stokes equations.

Boltzmann equation
A fluid can consider like a continuous medium when the fluid is dense, this meaning that the particles are behind the other particles and there isn't space between they, in this case the fundamental equations that conduct the fluid's evolution are of kind Navier-Stokes (NE).One form of determinate if the continuous medium is acceptable, is through of Knudsen's number ( ) λ Kn= l , which is defining like the relation between the free mean trajectory λ (mean distance that a molecule travels through before collisions with other molecule) and the characteristic length l; the model continuous medium is acceptable for a rank of ≤≤ 0.01 Kn 1 .In other words the Knudsen's number should be less that unit, in this form the continuous hypothesis will be valid.
In other way in the real fluid, the particles have strong and continuous interactions.In a gas named rarefied, the particles move during many time like free particles, except when collisions occurs of binary body.
When the Knudsen's number is in the rate of ≥ Kn 0.01 , the fluid could be considered like flow of rarefied gas.
The LBM has its basis in the approximation of a fluid like rarefied gas (or diluted gas) of particles.The rarefied gas can be describing by Boltzmann equation that was derivative the first time in 1872 by Ludwig Boltzmann (Duderstadt & Martin, 1979;Schwabl, 2002;Chapman & Cowling, 1970).The NS equations can derive of Boltzmann equation in a limit.Any solution about Boltzmann equation needs an expression for collision term Ω(f).The complexity of it, carry the search of simple models of collision processes, it will permit to make easy the mathematical analysis.Perhaps collision model more known was suggested simultaneously by Bhatnagar, Gross and Krook (Bhatnagar et al., 1954) and it is known like BGK:

The
Where υ is an adjusted parameter and ( ) eq fr , v , t   is the local thermodynamic equilibrium of distribution.This simplification is called the "single-time-relaxation" approximation, already that the absence of space dependence, is implicating that () () . This model keeps important properties of collision term in Boltzmann equation.For example, satisfies theorem H and obeys the laws of mass, momentum and energy conservation (for this reason, it keeps the structure of hydrodynamic equations), it keeps the collision invariants The collision term with the model BGK is writing the following form Without knowing the form of Ω(f) there are however several properties which can be deduced.If the collision is to conserve mass, momentum and energy it is required that () The terms i , i 0,..., 4 are frequently called the elementary collision invariants since Any linear combination of the ψ i terms is also a collision invariant.
The ordinary kinetics theory of neuter gas, the Boltzmann equation is considered with collision term for binary collisions and is despised the body's force F α .This simplified Boltzmann equations is an integro -differential non lineal equation, and its solution is very complicated for solve practical problems of fluids.However, Boltzmann equation is used in two important aspects of dynamic fluids.First the fundamental mechanic fluids equation of point of view microscopic can be derivate of Boltzmann equation.By a first approximation could obtain the Navier-Stokes equations starting from Boltzmann equation.The second the Boltzmann equation can bring information about transport coefficient, like viscosity, diffusion and thermal conductivity coefficients (Pai, 1981;Maxwell, 1997).

Equilibrium distribution function and theorem H
The equation (1) with the collision term described for binary collisions, it isn't lineal, is for that reason that the solution is very difficult.Nevertheless, exists a solution for Boltzmann equation, it isn't trivial and is very important and is known like distribution function Maxwellian.For this case the Boltzmann equation presents a non reversible behavior and distribution function lays to distribution Maxwellian, this represent the situation of an uniform gas in stationary state.
For derivate the distribution fuction Maxwellian, supposes the absence of external forces F α and uniform gas, the distribution function is ( ) ff v , t =  .Whereas these conditions, Boltzmann equation (1) with binary collision term is Where takes on that particle is spherical and symmetrical, for that reason the integral with regard to the angle ε can be evaluated immediately.
The function H is Equation ( 6) can differentiate with regard to time ( ) If considers inverse collisions, it means particles with speed v  and v'  collided and moved outside with speeds v  and v'  .The result of these collisions will be ( )  The equation ( 12) is Maxwell distribution function for a gas and describes the equilibrium state of distribution function f.The constants can find if replace equation ( 12) whit conservation laws obtains the common form of Maxwell's distribution function ( ) The theorem H shows that any distribution function f lays out its equilibrium state eq f for long times.
If a system is not uniform, it is not in thermodynamic equilibrium then can obey law of Maxwell's speeds distribution.However, if "equilibrium absence" is not big, can considered like good approximation, all little volume (microscopic scale) is in equilibrium (considered like subsystem).This is for two reasons.First little portions of gas contain a big number of molecules.Second the necessary time for established the equilibrium in a little volume is brief in comparison with necessary time for that transport processes get equilibrium in little volume with rest of system (it is true when concentration, temperature, etc. gradients are not too much big).In consequence, can suppose that is local thermodynamic equilibrium so speed distribution in any volume element (macroscopic) of medium is Maxwellian, although density, temperature and macroscopic velocity change the position (Duderstadt & Martin, 1979;Schwabl, 2002;Bhatnagar et al., 1954;Pai, 1981;Maxwell, 1997;Succi, 2001;Succi, 2002;Cercignani, 1975;Lebowitz & Montroll, 1983).

From the Boltzmann equation to the lattice Boltzmann equation
Lattice Boltzmann equation can be obtain through two ways, first is through of "cellular automaton" and second starting from Boltzmann equation, it was review previously, for carries out derivation of Boltzmann's lattice equation is necessary the space time discretization.Immediately presents brief description of second way, it shows by pace series.
The derivation begins when the Boltzmann equation with the model BGK is writing like ( ) In the equation ( 14), is distribution function of only particle, ξ  is microscopic velocity, λ is the relaxation time due to collision, and g is the Boltzmann-Maxwellian distribution function (fM), is important mention that collision term has been transforming in accordance with equation (2).

Hydrodynamic properties of fluid, density ρ, velocity u
 and temperature T, could be calculated starting from momentums function f.For quantify the fluid's temperature, uses the energy density ρє (He & Luo, 1997).

Time discretization
Equation ( 14) it formulates in form of ordinary differential equation (ODE) Where is substantial derivative length wise of microscopic velocity ξ


. Equation ( 18) that mentioned is a lineal ODE of first order, which can be integrated over a time step of δt.
Using the integral factor method is obtaining Assuming that δ t is small enough and g is smooth enough locally, and neglecting the terms of order of ( ) O δ or smaller in the Taylor expansion of the right hand side of the equation In the equation ( 20) t λ τ δ ≡ is the dimensionaless relaxation time; this equation is very similar to lattice Boltzmann equation.For obtaining the equation is necessary the space velocities discretization, too equilibrium function g it could be consistent with Navier-Stokes equations.The equation ( 20) is the evolution of distribution function f with discreet time (He & Luo, 1997;Maxwell, 1997).

Approximation of equilibrium distribution
A point of view very import, is the obtaining of equilibrium distribution function, it was mentioned before, of this depends obtaining appropriately the Navier-Stokes equations.Maxwell's distribution that is employing like equilibrium distribution g for unitary particle mass and "D" dimensions is It could be expressed the following form () ( ) The function "g" expands through of Taylor series in u  until third order and is necessary eliminate these For little velocities (or little Match numbers), this is an exact approximation.

Discretization of the velocities
For discretization of the velocities, it will be -∞ a +∞ in both directions "x" and "y" for specific case in two-dimensional model (D2Q9), it will unroll in this part of chapter.The particle momentums of distribution function are very important, because of this depends the consistence of (N-S) equations.in the same way, the isotropy is keeping during the discretization, it is an important property in the symmetry of NE equations, of this form, lattice will be invariant for problem rotations.
Is important mention that an isothermic model is only necessary first momentum.These can be describing by the equation ( 24) two dimensions like The integral of equation ( 25) has the next form ( )


is the momentum function, it contains powers of velocities components.For recuperate model D2Q9 for square lattice is using the system of Cartesian coordinates and ( ) in previous equation, x ξ y y ξ are components of "x" and "y" of velocity ξ  .The integration of equation ( 25) using these values gives us the following equation for the equilibrium distribution function for two dimensional, 9-velocity LBE model: The w α are 0 where c is the sound speed of the model (ThÄurey, 2003).

From lattice Boltzmann equation to the Navier-stokes equation
The popular problems of kinetics theory is the derivation of hydrodynamic equations, in certain conditions, solution of ( )   fr , v , t transport equation is similar the form that can relate directly to continuous or hydrodynamic description.In certain conditions the transport process is like hydrodynamic limit.In 1911 David Hilbert was who proposed the existence Boltzmann equations solutions (named normal solutions), and these are determinate by initial values of hydrodynamic variables it return to collision invariant (mass, momentum and kinetics energy), Sydney Chapman and David Enskog in 1917 were whose unrolled a systematic process for derivate the hydrodynamic equations (and their corrections of superior order) for these variables.
In spite of have been proposed many approximated solutions to Boltzmann equation (including the Grad's method of 13 moments, expansions of generalized polynomial, bimodal distributions functions), however the Chapman-Enskog is the most popular outline for generalize hydrodynamic equations starting from kinetics equations kind Boltzmann (James & William, 1979;Cercignani, 1988).

Chapman-Enskog expansion
For show that ERB can use for describing the fluid's behavior, NS equations are derivate by process are named Chapman-Enskog's expansion or multi-scale analysis.It depends of Knudsen's number it was mentioned at the first part of this chapter; it is the relation between the free mean trajectory and the characteristic length.
For derivation of NS's equation, the Boltzmann equation divides in different scales for space and time variables.It has his basis in the expansion of parameterε , it is necessary for using the Knudsen's number.In general the expansion is truncating after second order terms.The following representation is for temporal variables The advection and diffusion are considered similar in space scale "x 1 ".While, the representation of differential operator is similar For consistent expansion, is necessary the second order term in space.Momentum equations of f are expanded directly an addition the next form The dependence of time f is on account of variables ρ,  u and T.
The NS equation can be recovered starting from analysis of lattice Boltzmann equation Expanding Equation (32) in both space and time up to second order yields ) The three scales from O (ε 0 ) to O (ε 2 ) can be distinguished in Equation ( 33), and are handled separately.In the following, subsequent expansions of the conservation equations will be performed.Giving as result the continuity equation to ε 0 , ( ) To ε 1 , equation ( 33 (37) (He & Luo, 1997;ThÄurey, 2003).

The model
The Lattice Boltzmann Method (LBM), its simple form consist of discreet net (lattice), each place (node) is represented by unique distribution equation, which is defined by particle's velocity and is limited a discrete group of allowed velocities.During each discrete time step of the simulation, particles move, or hop, to the nearest lattice site along their direction of motion, where they "collide" with other particles that arrive at the same site.The outcome of the collision is determined by solving the kinetic (Boltzmann) equation for the new particledistribution function at that site and the particle distribution function is updated (Chen & Doolen, 1998;Wilke, 2003).Specifically, particle distribution function in each site () , it is defined like a probability of find a particle with direction velocity e i  .Each value of the index i specifies one of the allowed directions of motion (Chen et al., 1994;ThÄurey, 2003).
For this work we use D2Q9 and periodic boundary conditions in the inflow and outflow plane and non-slip boundary (bounce-back) conditions on the walls and the porous matrix.Bounce-back conditions were used whenever the fluid hit a node of the porous matrix.Our porous media is represented by blocks that are projections in the plane of actual three dimensional geometries Stability is improved by considering the porous matrix as made out of these blocks and makes the code less noisy as well.To initialize the lattice, a constant body force (F) is used and acts during the simulations, which physically corresponds to a constant pressure gradient.In this work we focus only on externally applied pressure; namely we deal with pressure-driven flows.
Fluid flow in porous media is modeled here by using a modification of the Lattice Bhatnagar-Gross-Krook (LBGK) technique (eq.32).
Fluid density and velocity are calculated as follows:  ue .The method used in this work involves a modification of the LBGK, where instead of constant viscosity, the effective viscosity is used directly as a rate-dependent relaxation time parameter avoiding the calculation of the matrix elements for ij Ω (Boek et al., 2003;Ahronov & Rothman, 1993;Rakotomala et al., 1996).τ Is the relaxation time provided by the kinematic fluid viscosity.
To test our algorithm, simulations were run for fluids in a rigid pipe with n = 0.33, 0.56, 1.0, 2.0 and k = 0.001, 0.005, 0.5, 10.0 respectively, in a 200 × 200 lattice.The steady-state velocity profiles were reached after 30,000 times step.Dimensionless velocity profiles from the LBM and analytical solutions are shown in Figure 1.Our numerical results and analytical solutions are within a 0.9% error (the maximum value being for n>1).
Next step in our model is the simulation of the fluid flow in porous materials (or packed beds).These materials have a portion of space occupied by heterogeneous or multiphase matter.In studying flows through porous media, we extend the LBM so as to study arbitrary and random approaches.Even at a first glance, many porous media, particularly the natural ones, show that the distribution of pores is random (Christakos & Hristopulos, 1997), so randomness is a really important factor to describe natural and man-made porous materials; they usually possess formidably complicated architecture.In modeling, reasonable idealizations have to be assumed (Dullien, 1992;Gibson and Ashby, 1988;Telega and Bielski, 2003).

Velocity profiles
Simple predictions concerning the flow of non-Newtonian fluids through porous media (Bird et al., 1960) have long been provided.Though Darcy's law for the flow of Newtonian fluids through a porous medium agrees with a model of parallel tubes, analogous behavior is not obvious for non-Newtonian fluids.This model does not take into account the tortuosity of the flow path among others; an additional nonlinearity in this case.As a first step, by assuming that a porous medium can be viewed as a collection of parallel pipes, the relation between flux and force for a non Newtonian fluid can be approached by where q is the volumetric rate of flow per unit area.C is a proportionality constant, function of the consistency index k, the porosity ε and the effective permeability K.For n = 1, www.intechopen.comequation ( 39) reduces to Darcy's law.In practically all cases, eq. ( 39) is valid as long as the generalized Reynolds number; does not exceed some value between 1 and 10.Here, D ρ is a length dimension scale proposed by Collins (Geankoplis, 2000).
Our LBM are such that all experiments performed in this work are in the range gen Re 1 < and can provide some output in the relative importance of the different parameters (for the fluid, and the geometrical structure) with statistical validation.We fully comply with the empirical macroscopic scale relation (39) obtained experimentally.
Two different types of porous media are studied here: a) Arbitrarily generated porous media.In this case, the authors propose a set of porous media totally constructed by arbitrary choice, fulfilling the following important conditions.At first percolation is guaranteed.Secondly, well defined, interconnected channels are constructed and finally low tortuosity is provided.b) Randomly generated porous media.In this case, the medium is generated by means of a random configuration.The Box-Muller method for generating standard Gaussian pseudorandom numbers is used to obtain the positions of the seeds in the solid matrix, 3 X 3 lattice sites nodes are defined as blocks around each seed; thus providing the porous matrix.This simplified geometry substitutes the projections on a plane of more complex and realistic cases (throats channels, chairs, cylinders, etc) improve stability and avoids noisy results.
The core of the Box-Muller method is a transformation that takes as inputs random variables from one distribution and as outputs random variables in a new distribution function.It allows us to transform uniformly distributed random variables to a new set of random variables with a Gaussian distribution.We start with two independent random numbers, x 1 andx 2 , which come from a uniform distribution (in the range from 0 to 1).Then, apply the to get two new independent random numbers that have a Gaussian distribution with zero mean and a standard deviation of one.This method of generating random porous media produces similar pore space characterization that several natural media.
Moreover, as it is well known for a Newtonian fluid that standard Darcy law can be used in order to obtain permeability.Namely ( ) Here K is permeability and µ is the viscosity for a Newtonian fluid.In the case of a non-Newtonian fluid with power law viscosity, we introduce the apparent viscosity (instead of the standard viscosity) which is given by the following equation.

Local effective permeability
We work with porosities 47.4%, 51.9%, 68.7% and 71.9%, randomly and arbitrarily generated, outputs of our experiments are presented (see figures 2-4).This is done for n=2, n=1, n=0.529 and different pressure gradients.Actually, if n=1, the permeability bands do not depend on the pressure gradients as expected in natural systems for porosities 50% or larger.In these Figures, the pore space characterization is shown at the upper left corner of the figure.Secondly, right below, the flow paths in each particular case are shown (from the numerical data, vuggy zones can be detected, among others).Finally, in the right part of each figure, the permeability bands produced by our LBM are introduced.Although qualitative in part, important contributions of our paper are the correlations of the band structures with different pressure gradients (F), fluid parameters (n, the index number), geometrical factors (porosity and tortuosity), and randomness.
In Figure 5, we explore the role of different porosities and pressure gradients in arbitrarily generated porous media, for a non-Newtonian fluid (n=0.529).The progressive introduction of high permeability bands as a natural consequence of the geometry, fluid, and flow parameters shows interesting patterns, as can be observed.
Uncertainty is well recognized as an important factor in natural porous media (in fact, it is well known that effective permeability is subject to greater uncertainty than porosity) and numerous studies have employed random methods to model Newtonian flow in subsurface porous media by assuming a given effective permeability probability density function.
It is important to construct models that are able to closely mimic the heterogeneity of actual porous media and sufficiently efficient to allow simulation of flow and transport phenomena.To predict the network flow at core scale (for instance, in hydrocarbon reservoirs, packed beds and aquifers), we propose to construct the permeability probability distribution within our model.Our results are shown in Table 2 in this work.Although it is widely understood that the selection of a particular effective permeability probability density will markedly influence simulation results in applications, only a few studies (this paper, among others) describe the manner in which to construct these effective permeability probability density functions from mesoscale information.
We correlate the effect of randomness, porosity, and fluid parameters with permeability fields and probability distributions predicted by our model, thus improving our understanding of heterogeneous media for applications to natural systems.Normally, at core scales for reservoirs and aquifers; among others, the unknown permeability distribution in the subsurface on all length scales is much needed for practical goals (Sitar et al., 1987;Cooke et al., 1995).
Normally, first step in estimating permeability from thin-sections of natural media is to convert representative digital images into binary images.These binary versions are used for porosity estimation and as conditional data for the stochastic pore-structure realizations. www.intechopen.com Simulation of Non-Newtonian Fluids Through Porous Media 89 Flow simulations are then conducted on each pore-structure realization.Our approach then provides the local effective permeability estimate for each bed type (see for instance Figures 2 to 5).We perform experiments for non-Newtonian fluids and obtain permeability bands (see Figures 2-5) for the following power-law fluids A) n = 0.529, B) n = 1 and C) n = 2.This is done in the following cases: a) arbitrarily generated porous media and b) randomly generated porous media, porosities 47.4%, 51.9%, 68.7% and 71.9%, and a number of pressure gradients, see Figures 2 to 5. Also arbitrarily generated porous media, flow paths, and the obtained permeability fields for a number of porosities, pressure gradients and n=0.529 are introduced in Figure 5.We then perform a null hypothesis analysis over the data by Kolmogorov-Smirnov contrast, obtaining the best fitting effective permeability distribution for each case, as shown in Table 2. Kolmogorov-Smirnov contrast (a nonparametric null hypothesis analysis) is a technique where two data distributions are compared (the one to be tested and another one hypothetically true) and they are accepted to be statistically the same, provided the maximum distance between both is below a certain threshold.

Conclusions
The 2-D permeability variations considered in this work are applicable to real contaminated aquifers where the length and width of high-permeability inclusions are large relative to their height (or thickness).For example, fluvial deposits as ''cut-and-fill'' (e.g., scour pool and dissection elements) and accretionary (e.g., horizontally bedded gravel sheets) elements can have high length-and width-to-thickness ratios.The 2-D permeability variations considered in this work are not applicable to more complex contaminated aquifers, e.g., where high-permeability inclusions do not have large length-and width-to-thickness ratios, or where tortuous 3-D flows cause streamlines to twist.
Effective permeability distributions have been used to study geological phenomena for petroleum technology, among others.In this paper, a study of how different geometrical features and rheological parameters affect local permeability of a porous medium has been conducted in order to enhance the knowledge of the local permeability distribution and its importance when developing global permeability models (Lundström et al.,2004;Lundström & Norlund, 2005).We conclude that randomness, porosity, and the power law index for non-Newtonian fluids play an important role in the local effective permeability distribution obtained for each experiment (Normal, Gamma, or any other).Our detailed results are presented for a range of porosity from 47.4% to 71.9% in this paper.We are able of using our approach for porosities as low as 10% , but then the average relative error on n (comparing the experimental with the value obtained by LBM) can reach a maximum of 6.8% (for n=0.529).See Table 2.
We consider several heterogeneous porous media representations to illustrate flow performance; all of them are two-dimensional.A number of high-permeability inclusions within low-permeability zones are obtained from a variety of fluid characterizations (see Figures 2-5).The effect of porosity, randomness, fluid and flow parameters are correlated with permeability bands in a number of cases; original outputs shown in this paper.
It is intriguing to see how Newtonian fluids always produce symmetrical bands and high permeability bands included into small permeability bands with several layers provided porosity is large.We can also remark that for n<1, bands of low permeability are more numerous than in the case of n>1.Some throats can be formed if we introduce randomness and high porosity, see Figure 4, in the case of n<1 and high pressure gradients.Additional anisotropy appears for small porosity, see Figures 2.
Streamlines determined by numerical simulation are obtained and shown in Figures 2-5.From these, and our numerical results (rough input data to construct Figures 2-5), it is clear that flow is mainly parallel and uniformly distributed in the low-permeability zone.However, in the region containing the high-permeability inclusion there is very little flow in the low-permeability zone and almost all flow is focused in the high-permeability inclusion.In other cases, a number of high-permeability inclusions within a low-permeability zone are produced.In regions between high-permeability inclusions, most flow is uniformly distributed, see Figure 5.In regions containing high-permeability inclusions, most flow is focused in these inclusions and is parallel to the boundaries.From our numerical results and the path flows, vuggy zones can also be easily identified.
The effective permeability data (input for our figures), represented by path flow, show that the "vuggy" features are always in the neighborhood of the walls distributed around a central "non-vuggy" zone.These results show particularly low effective permeability near the entrance and towards the end of the tunnel.It is clear that a zone of high effective permeability can be found towards the end of the tunnel, especially for Non-Newtonian fluids and low effective permeability near the tunnel walls.Similar studies for effective permeability maps, as functions of the distance from the tunnel wall and the distance from the entrance, for sandstone perforated with underbalance have already been reported (Karacan and Halleck, 2000).
It is straightforward to see, from our experiments, that if n (fluid parameter, power law index) increases, the number of homogeneous effective permeability zones decreases; namely the local probability data take values in a smaller set of numbers.Besides, for Newtonian flows, the effective permeability field is always symmetric for large porosity, no matter the detailed structure of the porous media.For non-Newtonian fluids, this symmetry may be affected by an angle sweeping the horizontal axis, provided the porous medium is not randomly generated (here, symmetry in the local probability field is lost) and most flow is globally oriented at a clear angle to the horizontal direction.
Pseudoplastic fluids in randomly generated porous media generate effective permeability fields, where zones of constant values are oddly shaped (bottle necks, among others, that can be considered scale up, equivalent to geometrical features, result of our mesoscale experiments), totally abnormal behavior can also be obtained for a different set of parameters.There seems to be a critical value for porosity to keep symmetric effective permeability fields.Our results can be used as models for materials with several zones where different fine erosion and deposition processes occur and provides an original approach to learn from the effect of the fluid parameters (n, the power law index) in the flow, provided the remaining features are constant.
Finally, in what follows, we learn how our LBM (mesoscale information) can be used, by the inter-relationship of geometrical factors, randomness and fluid parameters, to produce permeability (as probability laws) distributions in order to predict scaled up subsurface permeability bands.
In our experiments, normal distribution for effective permeability is obtained mostly for randomly generated porous media provided the fluid is Newtonian or shear thinning, see Table 1.Eventhough this distribution is not expected for arbitrary generated porous media, it may appear provided the fluid is Newtonian.It should be remembered that the use of normal distribution can be theoretically justified in situations where a large number of effects act additively and independently together.The porous media, LBM and fluid factors used in these experiments can be a first approach to the effective permeability distribution of caffeine and testosterone solutions in silicone membrane (Khan et al., 2005).
Provided porosity increases, our results show that for Newtonian and shear thickening fluids, regardless the way the porous matrix (randomly or not) was generated, the effective permeability distribution tends to be Weibull.It is possible to conclude that the experiments shown in Table 1 as Weibull effective permeability distributions are describing systems involving a weakest link and/or non-linear effects derived from the fluid viscosity.Statistics of the mechanical and failure properties on the grain scale are often assumed to follow the Weibull distribution; here, the significant influence of microcrack length statistics has been emphasized (Wong et al., 2006).A Weibull distribution can accurately describe experimentally determined time trends of the infiltration rate (Faybishenko et al., 2003).
Infiltration experiments conducted on packs of rocks show that fluxes may stabilize into an exponential distribution (Tokunaga et al., 2005).Exponential distribution for effective permeability has also been reported in the analysis of heterogeneities in porous media (Savioli et al., 1996).In our simulation, the effective permeability distribution is exponential only for arbitrarily generated porous media and shear thinning fluids, no matter the porosity used in the simulation, see Table 1.They are fit to be used as models for petroleum processes.We believe that arbitrarily generated porous media may provide the deterministic background to obtain exponential effective permeability distributions since these are found in situations where certain events occur with a constant probability per unit distance.Clearly, for non-Newtonian fluids, non-linear effects are introduced; these are small if n tends to one.
The gamma distribution has been used as a statistical representation of surface roughness for flow in unsaturated fractured porous media and tortuosity distributions in porous media (Or and Tuller, 2000;Lindquist et al., 1996).In our simulation, the effective permeability distribution is gamma for shear thickening fluids and rather low porosity, regardless the way the porous matrix was produced, see Table 1.So the gamma effective permeability distributions found in our experiments can be considered applicable to particle tracking.This process has been modeled as the sum of elementary steps with independent random variables in the sand.
In our experiments, the effective permeability distribution is Generalized Extreme Value distribution for arbitrarily generated porous media, porosity 71.87%, n=2.This seems to be a very particular case, probably coming from a strong non-linearity.This is the limited distribution of the maxima of a sequence of independent and identically distributed random variables.Extreme hydrometereological phenomena follow an extreme distribution (Gutiérrez-López et al., 2005).At constant potential, the maximum peak currents in different time intervals during a potensiostatic test follow extreme distribution (Zuo et al., 2000).
Summarizing, this paper has presented a method to study scale up processes (mesoscale to core scale) of flow in porous media from geometry, fluid parameters, stochastic realization and LBM; namely, parameters such as tortuosity (arbitrarily generated porous media always have a clearly smaller tortuosity than randomly generated porous media), randomness in void space, geometrical distribution, porosity and consistency index of the non-Newtonian fluid, among others.This model is ready to predict as an output from thin small samples, effective permeability bands and the corresponding subsurface probability distributions scaled up for subsurface reservoirs, aquifers, etc.

Acknowledgment
The first author wishes to thank CONACYT and UNAM for financial support.
, v during the time t.The micro state is given by the complete set { } αα   ) Eule'r equation.Where

Fig. 1 .
Fig. 1.LBGK velocity profiles and Analytical Navier Stokes velocity profiles in a rigid pipe.
permeability for non-Newtonian fluids (with power-law viscosity) is calculated in this work as follows
Time t represents local relaxation that is very quickly in fluid for collisions.The diffusion processes are time scale t 1 .In this only space expansion is considerate, it is represented in the next first order expansion.

Table 1 .
Best fitting permeability probability distributions for number of experiments. www.intechopen.com

Table 2 .
Average relative error as a porosity function and the experimental value of n in our method.