3D lattice Boltzmann simulation of fluid flow in fractured rock

3D lattice Boltzmann simulation of fluid flow in fractured rockT.Miyoshi,T.Matsuoka,S.Murata&Y.AshidaKyoto University,Kyoto,JapanABSTRACT:It is important to estimate fluid flow in reservoir rocks for an enhanced oil recovery and evaluation of the productivities of the fracture reservoir.In this study,the Lattice Boltzmann Method(LBM),which is a new numerical computation method for the fluid flow,was applied to fracture flow.Consequently,we concluded that the flow obeys the cubic law even in the case of a very complex surface when the parallel fracture surfaces are the same shapes.The deviation from the cubic flow,however,is controlled by the standard deviation of the aperture of the fracture.This means that the variance of the topology of the fracture surface is not an important parameter compared with the variance of the aperture width.1INTRODUCTIONWe have conventionally used Reynolds equation for evaluating fluids flow behavior in single fracture,but we cannot evaluate the exact channel flow.In a simple approach using Reynolds equation,the vertical flow in the fracture cannot be evaluated.To overcome this problem,the Lattice Boltzmann Method(LBM)(Qian et al.1992,Rothman&Zaleski1997,Succi2001, Wolf-Gladrow2000,Y oshino2000)has been applied in this study.LBM has an advantage that we can evalu-ate the3D flow behavior with the complex boundary conditions.In this study,we developed2D&3D LBM programs and applied them to the fracture generated by the fractal algorithm.We have done the following simulation studies.The fluid flow generally becomes less than the theoretical quantity when the fracture aperture is narrower.This phenomenon is shown in this simulation and the deviation from the cubic law is discussed in detail.The two parameters of fracuture topology influence this deviation.2BASIC THEORIES FOR LBM2.1Lattice Boltzmann methodA finite difference method or a finite element method have conventionally been used for the simulation of fluid behavior.These methods are based on Navier-Stokes equation.In contrast to these method,the LBM approximates fluid as particles.The collision and movement of particles and boundary conditions are expressed by density distribution functions.And the fluid velocity,density and pressure(macroscopic parameters)are calculated from the distribution func-tions(microscopic parameter).The distribution func-tions must be set as if the collision,movement,and boundary conditions retained the character of con-tinuation of fluid.Various models of different lat-tice or collision rules are suggested,but we adopted the BGK(Bhatnagar-Gross-Krook)model(Bhatnagar et al.1954)and set9velocity model for2D and15 velocity model for3D.2.2Basic equationsThe concept of the LBM in the15velocity model for 3D is as follows.At first the three dimensional space is divided into unit lattices.It is assumed that the particles are on the lattices,and the number of particles cor-responding to a particle density are expressed by the distribution functions.Furthermore,the virtual par-ticles on eachlattices have15velocity,as expressed in Equation1.Figure1shows the velocity of each particle.In addition,all particles exist on the lattices and there are no particles on the way of between lattices in each calculation time step.The distribution functions according to each15velocity are given on all lattices. It is very important for the distribution functions to change with time steps in the LBM.Figure 1.The three dimension lattice is shown.The virtual particles move to each 15direction.There are four parts in this method.The first step is the collision part.The evoluation of the distribution function f i (x ,t )with the velocity c i at the point x and time t is calculated by the followingequationswhere x is a lattice spacing, t is a time step,τis a single relaxation time and f eq is an equilibrium distribution function that is given from fluid macro-scopic velocity u and density ρusing Equation 3.The equilibrium disutribution function is a particle density distribution that reached equilibrium in the finite space likelattices.The second step is the movement part.The particles on each lattice move in 15directions according to par-ticle velocity c i when a time step advances.The number of particles which move in each direction is defined as the distribution function.The third step is the evalu-ation of the boundary conditions.We will describe this in detail in the following section.The last step is the macroscopic part.The fluid velocity,density and pressure,which are macroscopic parameters on each lattice,are calculated from the distribution functions (Eq 4to Eq 6).These distribution functions are given from the movement step.When these four steps arerepeated,the evoluation of fluid flow can becalculated.2.3Boundary conditionsAll the distribution functions are not known on the sur-face of a boundary wall because the particles does not move from outside of the fluid area.So the unknown distribution functions must be determined from the known distribution functions.Several boundary con-ditions are proposed.We use the Inamuro ’s method in this study (Inamuro et al.1995).When a lower half of Figure 1is a wall,the distribution function f 3,f 8,f 9,f 11and f 14having speed of c y >0are unknown.Then each distribution function is guessed as Equation7.where u w ,v w and w w are wall velocity to each direc-tion.The wall velocity is zero in thissimulation.In Equation 7the unknown parameters,ρ ,u ,and w are given as follows,since the fluid velocity must be equal to the wall velocity and three equations on vel-ocity are made according to each direction.Moreover,the fluid density on the boundary lattices ρw is unknown quantity and is calculated by Equation 4.Equations 8–11are provided as a solution of theseequations.Next,for the boundary condition with an intersecting borderline and corner of a wall we assume a virtual boundary surface including boundary node and apply the above boundary condition on the virtual surface.A periodic boundary condition is adopted at the inlet and outlet surfaces.This boundary condition is set according to the pressure difference between the inlet and outlet.Supposing that the pressure differ-ence is generated to the x axis and the surface x =0is the inlet and x =L x is the outlet,the disutribution function f 2,f 8,f 10,f 11and f 13are unknown on the node of the inlet surface and f 5,f 9,f 12,f 14and f 15on the outlet surface.Then as f 2is known at the outlet surface,this function is introduced into the inlet side.That is,adding the pressure difference K ,f 2at the inlet side is given by Equation 12.In addition,the other un-known distribution functions obey Equation 13.The unknown functions at the outlet side are determined,as well as the inlet side (Eq 14&15).The pressure difference K is given by Equation16.3THE FLUID FLOW IN THE FRACTURE 3.1The deviation from the cubic law of fracturepermeability It is known that the flux through a single fracture is proportional to the cubic of mean aperture if a model is a parallel flat wall (Brown 1989).But as the aper-ture narrows,the flux does not obey the cubic law.We examined this result using the LBM.Since the fracture surface is generally fractal,a fracturesurfaceFigure 2.2D fracture model.This model is divided into 1000×100lattices.Figure 3.The deviation from the cubic law for the 2D fracture model.The deviation becomes larger according to narrowing the aperture.model was made by the fractal algorithm,as shown in Figure 2.The two fracture surfaces were parallel each other.The boundary conditions was set as follows.The pressure difference at two boundaries each other was set and the fluid flow in this direction was calculated.Next,the fracture aperture was changed by moving the upper surface upwards or downwards.The mean aperture is defined by Equation17.where d (x ,y )is the distribution of the fracture aper-ture,L x and L y are the length of a fracture to each direction.Now,when the flux Q is calculated by the LBM,the hydraulic aperture that conforms to the cubic law is determined by Equation18.where µis the viscosity of the fluid and p is the pressure difference to the x direction.If the ratio of d h /<d >is equal to one,this fracture model perfectly conforms to the cubic law.The cubic of this ratio means a deviation from the cubic law of fracture permeability.Figure 3shows the result of this simulation.In this Figure,the vertical axis is the cubic of the ratio (d h /<d >)3and the horizontal axis is the mean aperture normalized by standard deviation of fracture surface.It can be seen that the (d h /<d >)3deviates from the cubic law with the decrease in the aperture.This result is almost the same as Figure 10in Brown (1989).Figure4.The schematic drawings ofσf(a)andσd(b).σf is standard deviation of fracture height andσd is that of the aperture.3.2Parameters of the deviationIt is shown in the subsection3.1that the narrower the mean aperture,the larger the deviation from the cubic law.But if the model is a parallel flat wall,the devia-tion from the cubic law cannot be found theoretically no matter how narrow the aperture is.Then a roughness of fracture is considered as the parameter of the devi-ation.In this study,the standard deviation of fracture heightσf and the standard deviation of the apertureσd are used for the roughness parameter.Figure4shows schematic drawings ofσf andσd.The fracture model made by the fractal theory has almost the same value ofσf for the upper and lower surface.If the upper and lower wall are the same,σd is equal to zero because the aperture is uniform.3.3Analysis of2D modelsWe analyzed the deviation with the following2D models.Model2D-A:the upper and lower surface whose shape is a sine curve have the unifrom aperture (Figure5(a))Model2D-B:the upper surface is a sine curve and the lower is flat(Figure5(b))Model2D-C:the upper and lower surface are the same surface which is made by the fractal theory, and then,the upper wall is moved to the fracture length direction(Figure5(c))Several models of differentσf are made according to the amplitude of a sine curve in model2D-A.The bigger the amplitude,the larger the standard deviation of the surfaceσf.However,σd is equal to zero because the aperture is uniform in each fracture.In Figure6, the deviation of the models for differentσf is ploted as a function of the aperture.The legend“σ”areσf. The solid linewhose legend is“fracture”corresponds to that shown in Figure3.This figure shows that the deviation from the cubic law as a result of narrowing the aperture cannot be found even ifσf changes that.Figure5.Models for analysis of the parameter contribut-ing the deviation from the cubic law.The size of models is;333×33lattices(2D-A and2D-B);1000×100lattices (2D-C).Figure6.The lines with symbols are the rusult of model 2D-A.The solid line corresponds to Figure3.The legend “σ”is allσf.σd is zero.The deviation from the cubic law cannot be seen.In the model2D-B,the aperture changes along the profile because the lower wall is flat.Thenσf is equal toσd.The rusult of model2D-B is the three dashed lines with filled symbols in Figure7.Due not only to narrowing the aperture but also increasingσ,the deviation from the cubic law is increasing.In the model2D-C,at first the upper wall and the lower wall which are made by the fractal theory are the same.Next the upper wall is moved to the fracture length direction.Thenσf is fixed at12.6.The rough-ness of the aperture distribution increases according to the shift length of the upper wall.The result of2D-C is the three lines with unfilled symbols in Figure7.The deviation becomes larger with increasingσd.These things show that the standard deviation of the aperture distributionσd controls the deviation from the cubic law on the flux in the fracture.3.4ConclusionsThe flux in a fracture is proportional to the cubic of the mean aperture only when the aperture is uniform. That is,as the upper and lower walls are the same,the deviation from the cubic law does not occur even if the surface is too complex,but the deviation becomes larger withσd.The separation is strongly controled not by the roughness of the fracture surface but by theFigure 7.The result of model 2D-B &model 2D-C.The legend “σ”is σd .The dashed lines with filled symbols are the result of model 2D-B.The dashdotted lines with unfilledsymbols are the result of model 2D-C.The solid line corres-ponds to Figure 3.The deviation from the cubic law becomes larger with increasing of the σd .roughness of the aperture.As a conclusion,σd is a more important parameter to be considered in the deviation from the cubic law,and the deviation is larger with increasing roughness of distribution of aperture.3.53D fracture modelWe described the 2D models in the previous section but when we compare the results for this simulation to the experimental results,3D analysis is more valuable.In this section,we analyzed 3D fracture models using the 3D LBM,and compared them to 2D.In the 2D analysis,the 2D models were made cut-ting out from the 3D fractured model whereas in the 3D analysis the 3D fracture model is used as it is.We made various models for the same fractal dimension 1.2.And we paid attention to the σf and σd and exam-ined them as well as in the 2D analysis.Figure 8is a lower surface which is divided into 256×256×15lattices.And in Figure 9,the visualized 3D fluid flow is shown.The Axis value is a lattice number.The solid lines are streamlines and the gray plate is a section of the fracture.The vertical fluid flow which cannot be seen by Reynolds equation can be seen.3.6Analysis of 3D modelsTwo kinds of models are used for the 3D analysis.Figure 10shows the 2D section of the following models.Model 3D-A:(the tensile type fracture)The rough-ness of aperture σd is 0.65for all but the roughness of the fracture surface,σf ,is varied in each model.Model 3D-B:(the strike-slip type fracture)Using one of the 3D-A models,the upper surface is moved to the horizontal direction.The σd increases from 0.65according to the shift length but σf is constant.Figure 8.The lower surface of the 3D fracture model made by the fractal algorithm,which is divided into 256×256×15.Figure 9.The visualized fluid flow through the 3D fracture.A solid line is a streamline for fluid flow.A gray plate is a section of the fracture.Figure 10.The schematic drawings of the 2D section for two fracture types.In model 3D-A,we made three models whose σf is 2.97,1.88and 2.34,respectively.And the standard deviation of aperture σd is about 0.65.(d h /<d >)3for these three models are calculated using the 3D LBMFigure11.The model3D-A.The legend“σ”isσf.Theσd is constant.The trend of the deviation is the same in eachmodel.Figure12.The model3D-B.The legend“σ”isσd.Theσf is constant.The deviation becomes larger according to theσd.and plotted as a function of the aperture normalized by the standard deviation of the aperture(Figure11). Even if the roughness of surface(σf)is different, the trend of the deviation from the cubic law is almost the same.This result is almost the same as the2D analysis.In model3D-B,the surface which is one of model 3D-A(σf=1.88)is not changed as the shape,but it is moved in the direction of fluid flow.This model is assumed to be the fracture with strike-slip.Firstσd is 0.65,and it increases according to the shift length.The lower surface is fixed,and the upper surface is moved with X to the x direction.Then the upper wall does not exist at0≤x≤X.So using the periodic boundary condition,the empty area is filled by the upper surface data at L x−X≤x≤L x.In this condition,σd is maximum as L x/2shifts because the two surfaces do not fit.In this simulation,the shift length is4–8%of the fracture length.In the Figure12,the legendσis allσd,andσf is equal to1.88.As the fracture surface is almost equal,the deviation of the cubic law is larger with the increase inσd.3.7ConclusionsThis is the same result as the2D analysis.The deviation from the cubic law is almost the same whenσd is fixed andσf is changed.Whenσf is fixed andσd is changed, the separation becomes larger asσd increases.4CONCLUSIONSIn this study,we confirm that the lattice Boltzmann metheod is valuable for the simulation of fluid flow through the reservoir rock like the porous media or ing this method,the fluid flow in a complete3D space can be calculated more correctly.From our study,we concluded that when the paral-lel fracture surfaces are the same topology,fluid flow obeys the cubic law even in the case of a very complex surface.This means that the variance of the topology of the fracture surface is not an important parameter, compared to the variance of the aperture.The deviation from the cubic law is controlled by the standard devi-ation of the aperture of the fracture.In the tensile type fracture,the permeability will be near to that estimated from the cubic law because the upper surface and lower surface fit each other very well whereas in the strike-slip type fracture,the permeability will become less because the distribution of the aperture is disturbed. Therefore at the tensile type fracture,more amount of fluid can flow than the strike-slip type fracture which has the same mean aperture.REFERENCESBhatnagar,P.L.,Gross,E.P.&Krook,M.1954.A model for collision processes in gases.I.Small amplitude processes in charged and neutral one-component systems.Phys.Rev.94:511–525.Brown,S.R.1989.Transport of fluid and electric current through a single fracture.J.Geophys.Res.94:9429. Inamuro,T.,Y oshino,M.&Ogino,F.1995.A non-slip bound-ary condition for lattice boltzmann simulations.Phys.Fluids7:2928–2930;Erratum:8(1996):1124.Qian,Y.H.,d’Humieres,D.&Lallemand,ttice BGK models for Navier-Stokes equation.Europhys.Lett.17:479.Rothman,D.H.&Zaleski,ttice-gas cellular automata-simple models of complex hydrodynamics.Cambridge:Cambridge University Press.Succi,S.2001.The Lattice Boltzmann Equation for Fluid Dynamics and Beyond.New Y ork:Oxford University Press.ttice-gas cellular automata and lat-tice Boltzmann models.LNM1725Berlin:Springer.Y oshino,M.2000.Numerical analysis of transport phenom-ena in porous structure by the lattice Boltzmann method: Ph.D.thesis.Kyoto University.。