- Email: [email protected]
Simulation of uncompressible fluid flow through a porous media
Chaos, Solitons and Fractals 39 (2009) 1753–1763 www.elsevier.com/locate/chaos
Simulation of uncompressible fluid flow through a porous media A. Ramı´rez
a,* ,
J.L. Gonza´lez a, F. Carrillo b, S. Lopez
c
a
b
Instituto Polite´cnico Nacional (SEPI-ESIQIE-IPN), Unidad Profesional Zacatenco, Laboratorio de Ana´lisis Met. (Edif. ‘‘Z’’ y Edif. ‘‘6’’ P.B.), Me´xico City, Mexico Instituto Polite´cnico Nacional (SEPI-CICATA-IPN), Unidad Altamira Tamaulipas, Me´xico, Mexico c Instituto Mexicano del Petro´leo (I.M.P.-D.F.), Me´xico, Mexico Accepted 18 June 2007
Abstract Recently, a great interest has been focused for investigations about transport phenomena in disordered systems. One of the most treated topics is fluid flow through anisotropic materials due to the importance in many industrial processes like fluid flow in filters, membranes, walls, oil reservoirs, etc. In this work is described the formulation of a 2D mathematical model to simulate the fluid flow behavior through a porous media (PM) based on the solution of the continuity equation as a function of the Darcy’s law for a percolation system; which was reproduced using computational techniques reproduced using a random distribution of the porous media properties (porosity, permeability and saturation). The model displays the filling of a partially saturated porous media with a new injected fluid showing the non-defined advance front and dispersion of fluids phenomena. Ó 2009 Published by Elsevier Ltd.
1. Introduction The analysis of the fluid flow behavior in a porous medium involves the solution of some mathematical equations that must be calculated simultaneously to solve transport conditions for an unsteady state. For engineering problems, the solution is quickly found when the problem is related for a homogeneous medium and the boundary conditions are easy to define and represent. The equation of continuity for fluids can be solved using a finite difference or finite element method on rectangular grids, whose solution is easy to program using nested loops algorithms. But the problem becomes complex when a porous medium is analyzed. Because the properties involved in the transport of fluid flow are not constant. Averages of porosity, permeability and saturation conditions are frequently taken to represent the behavior of a system of porous media especially for analysis on macro or mega-scales; nevertheless, these methods not always represent in detail the behavior of the fluid flow inside. For example, in natural oil reservoirs with heterogeneous properties, there are irregular porous regions bordered by regions with minor porosities or permeabilities, for this reason they exist in the world as natural confined volumes. They can be formed by groups of big rocks with *
Corresponding author. Tel.: +55 57 29 6000x55127. E-mail addresses: [email protected], [email protected] (A. Ramı´rez).
0960-0779/$ - see front matter Ó 2009 Published by Elsevier Ltd. doi:10.1016/j.chaos.2007.06.105
1754
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
available space between them or long granular regions that allow the residence of a fluid inside under high pressure conditions. Here, one of the biggest problems is to describe the geometry and properties of these reservoirs. The shape of a rock or little grains and terrain geometry, all of them are irregular and cannot be analyzed or described using simple Euclidean or analytical geometry. This means that they cannot be represented as spheres, squares, ellipsoid or any regular geometrical form. Furthermore, in the real world, natural reservoirs are absolutely saturated and due to their irregular geographical position there is an irregular pressure distribution; this pressure is very high and produces oil inside springs up to the surface. A very popular kind of models to simulate fluid flow in anisotropic materials [14,5,11] was developed at the beginning of the 1940s with fractals and stochastic theories and it is called the ‘‘random walker’’ or ‘‘drunken walker’’. In this fluid flow through the porous medium is simulated considering a relative probability for following a non-deterministic path. Some of these models have been used to represent some natural phenomenon like the cracking propagation, the randomly way of rays, etc. [4,5]. The first step to represent fluid flow systems is to define the available spaces of the medium, the outlets, the inlets and their injection rates, so that the fluid injected into the system is named the new fluid or the invasion fluid by some authors [15,16] and the resident fluid is called the old fluid or the defender fluid. The new fluid begins to conquer the nearest points of the mesh and push out the resident fluid. These models were the platform to develop the percolation theory and are based on the employment of Monte Carlo methods to predict the probable fluid path [2,3]. These models have been studied mainly by oil companies to predict the oil behavior and evaluate and improve the efficiency of their recovery techniques. Another way to simulate fluid flow is to create a computational representation of the porous media and then calculate the fluid displacement by solving fluid flow equations as a function of the medium properties. This process is illustrated in this work. Here, a section of the porous media can be reproduced taken as a base the methods of Monte Carlo and random numbers generators included in logical algorithms in order to build a group of regions with an interconnected porosity in an original non-homogeneous medium, and these regions represent the possible channels where the fluid flow will be driven. These models as a function of the scale used represent the tortuosity of the paths and the hydraulic conductivity as a parameter to calculate fluid mobility (k), permeability of the medium, etc [6,7]. Random walk models can also be programmed to display a group of cracks or a group of rocks with interfacial spaces. The points where the cracks are joined or separated in many branches will be taken as the fluid way. The same consideration can be taken for the points that define the rocks borders; here the equation of conservation is used to solve the distribution of the fluid as a function of the new fluid received and then the excess of fluid on saturated capacity will be redistributed [10–13].
2. Mathematical model for PM Here a cylindrical core of sand berea is taken as a sample as shown in Fig. 1 with closed faces. A squared mesh of a finite difference method is used to represent a split of the porous media. The dimensions and volume that represent every node are the same and remain constant during simulation. This process is repeated for assigning the three most important properties of the porous media to every node (porosity, permeability and saturation). The 2D mathematical model used to mesh the sample is described in detail in a previous work [20]. The algorithms include a subroutine for
Fig. 1. Volume of the cores sample representing the mathematical model.
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1755
Table 1 Probabilities for defining porosity of the sample Class
Position I=1
I = nx
J=1
J = ny
1 2 3 4 5
40 20 15 15 10
40 20 15 15 10
20 20 20 20 20
20 20 20 20 20
/effe 4.00 2.00 6.00 8.00 0.00
Table 2 Probabilities for defining permeability of the sample Class
Position
1 2
K
I=1
I = nx
J=1
J = ny
50 50
50 50
50 50
50 50
5 2
reading and calculating the distribution of properties. The information for the distribution is shown in Tables 1 and 2 and was taken from real data of sandstone berea cores characterization. A reproduction of the core sample is generated using an algorithm based on a stochastic assignation value subroutine. Here every value of porosity is classified in classes (c) as is done for populations in a statistical analysis, according to their presence frequency in the real core sample. Then a probability to be present in PM is defined for the perpendicular faces. So that, the appearing probability of the porosity (PcI,J) is calculated using Eq. (1) as a function of the element position (I,J) in the sample. The process is illustrated in Fig. 2 and it is repeated for all the defined classes; the probability (Pc) is obtained by the sum of the particular probabilities due to this (PcI,J = 100). Instead the probabilities are compared with a generated random number (Z), and if the value of (Z) is found in the range (Pc 1 < Z < Pc) the value of the class (property) is assigned to the node I ðP /I¼1 þ P/I¼nx Þ þ nyJ ðP /J ¼1 þ P /J¼nx Þ nx : ð1Þ PcI;J ¼ 2 Some authors [1,9,17,18] have assumed empty core samples in their models according to their particular cases to be solved. Here, we consider a partially saturated medium with a non-equal probability for distribution of the saturation classes along the longitudinal axis ‘‘x’’ but constant for the transversal axis ‘‘y’’ in contrast with the profiles obtained for porosity and permeability that remain constant along both axes as is shown, respectively, in Fig. 3a–c. A general saturation profile as a function of the maximum pore volume defined (Vmax pore) is shown in Fig. 4a. Eq. (2) is used to compare the saturation values in order to display a new saturation profile as a function of the maximum pore volume defined. Additionally, isolated or non-accessible regions of the sample can be represented with blocked nodes. These isolated elements can be easily identified in a new equivalent representation as is shown in Fig. 4b. Isolated elements are commonly used to watch all the interconnected pores (effective porosity), V sat tot I;J ¼
V sat1;J : V max pore
ð2Þ
This process has a mathematical base that is described next. It is well known that fluid flow though porous medium is related with medium spatial variation; in consequence its behavior (displacement and advance front) cannot be accurately predicted using only deterministic methods. This work is focused on describing fluid flow in a complex porous media that produce non-defined advanced fronts as a consequence of a heterogeneous distribution of the properties involved. Uncertainties in the values of the porous media properties can be reduced but it can never be entirely eliminated; the theory of stochastic process provides a natural method for evaluating them, and a tool to simulate these cases using computational techniques. In the present work the properties of the porous media were treated as random space functions (RSFs), whose spatial correlations can be quantified by probabilistic distributions as was mentioned [20]. A better explanation is given next.
1756
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
Fig. 2. Flow chart of the subroutine to define PM properties.
Fig. 3. Profiles and distribution curves for core sample properties.
A random variable ‘‘U’’ can be defined by a set ‘‘X’’ of possible values (possible values of the properties) and a function of probability distribution designed by ‘‘P’’, which is defined on a collection of subsets of ‘‘X’’ denoted by ‘‘S’’; which means that it must be closed under countable unions and finite intersections, due to this ‘‘P’’ must satisfy certain
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1757
Fig. 4. Profiles (a) saturation and (b) isolated elements of PM.
conditions. For example if ‘‘P(X) = 0’’ this denotes an empty set, if (S1 S2) then P(S1) 6 P(S2) this implicates that the event ‘‘U 2 S’’. The probability that ‘‘U’’ lies between two different values that limits a range ‘‘u’’ and ‘‘v’’ can be written as follows: ‘‘P(u 6 x 6 v) = F(v) 6 F(u)’’. Here ‘‘F’’ is the cumulative density function; now if the difference between ‘‘u’’ and ‘‘v’’ is infinitesimal ‘‘du’’ the probability that ‘‘U’’ lies between ‘‘u’’ and ‘‘du’’ is ‘‘F(u + du) F(u) = dF(u)’’ and its probability of density function will be P(u) = dF(u)/du. Finally, the probability that ‘‘U’’ lies between an interval ½a; b can be written in an integral form as is shown: P ða 6 U 6 bÞ ¼
Z
b
P ðuÞ du ¼ F ðbÞ F ðaÞ:
ð3Þ
a
A procedure based on the comparison of the probabilities to lie between a pair of values was used to define and simulated the porous media properties in a simple computational algorithm. The node definition has two main purposes. The first is to store the particular properties of everyone and the second is to establish the starting conditions for the simulation of fluid flow. In other words it must be defined which nodes will work as inlets, pores and outlets and what are their corresponding injection rates. It means that it must be defined what nodes will participate in the fluid flow and what will be their roles during distribution of exceed of fluid.
3. Characterization of PM It is important to mention that there are two kinds of porosities. The absolute porosity (/abs) is the total empty volume of a material and the effective porosity (/effe), which is formed by all connected pore spaces from an initial point to an ending point in a PM. Both properties were characterized in this work and the procedures are the following. For the (/abs), samples of sandstone berea were cut, cleaned and polished resulting in thin plates which were prepared by introducing a colored resin into their porous structure, previously evacuated, and the porous regions are identified as the geometrical zones where the resin color remains. Thus, only the opened porous spaces are accessible to fluid. After this the samples were re-polished; then a set of five photographs of one face of five samples were snapped using optical microscopy. Finally, the relation between colored area with resin (Aresin) and the area of the photograph (Atotal) represents the measure of the (/abs) and it was counted using a software image analyzer (see Table 3). For the (/effe), a procedure according to the Archimides immersion method established on the recommended practices for core analysis [19] shown in Fig. 5 was followed to analyze and classify the samples. In consequence of the calculations it is possible to confirm (/abs > /effe), because there is a volume of unconnected spaces that contributes to the (/abs) but not to the (/effe). The results of the analysis are shown in Table 4. The average was (/effe = 4.10). This information is necessary to define the porosity data used to generate the simulation shown in Fig. 3a and eventually to solve the fluid flow equations.
Table 3 Probabilities for defining saturation of the sample Class
1 2
Position
S
I=1
I = nx
J=1
J = ny
80 20
20 80
80 20
20 80
80 90
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1758
Fig. 5. Archimides mercury immersion apparatus.
Table 4 Presence of (/effe) and (/abs) inside PM Class
(/abs)
(%)
(/effe)
(%)
1 2 3 4 5
10 12 14 16 18
10.00 26.67 33.33 20.00 10.00
0 2 4 6 8
10.00 20.00 40.00 15.00 15.00
4. Mathematical model for fluid flow displacement Here we used the finite difference method to solve Eqs. (4)–(7) for every node used in the discretization process at every iteration. The left terms of these equations represent the fluid flow and movement of the new fluid in a new node and the right terms the fluid flow and movement of the resident fluid ejected from the node (when their pore volume is absolutely saturated). Eqs. (4) and (5) are the continuity equations that can be solved as a function of the quantity of movement Eqs. (6) and (7) o qu dy dz ¼ qu þ ðquÞ dx dy dz; ð4Þ ox o ð5Þ qv dx dz ¼ qv þ ðqvÞ dy dx dz; oy o quV dy dz ¼ quV þ ðquV Þ dx dy dz; ð6Þ ox o qvV dx dz ¼ qvV þ ðqvV Þ dy dx dz: ð7Þ oy Continuity equations can be written as Eq. (8) and solved as a function of Darcy’s law (9). Here ‘‘k’’ is the conductivity that can be established as a function of the permeability and subjected to the conditions shown next. r qðxÞ ¼ gðxÞ; qi ðxÞ ¼ k s
ohðxÞ ; oxi
hðxÞ ¼ H B ðxÞ ðaÞ x 2 CD ; qðxÞ nðxÞ ¼ QðxÞ ðbÞ x 2 CN :
ð8Þ ð9Þ
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1759
These equations can be solved using a deterministic set for dependent variables ‘‘h(x)’’ and ‘‘q(x)’’ if ‘‘k’’ is known at every point of a grid, where the notations used means: x is the vector of space coordinates q is the flux discharged at point (x) $ is the grad operator with respect to (x) h(x) is the the hydraulic head kS(x) is the saturated hydraulic conductivity n(x) = (n1, . . ., nd)T is an outward unit vector normal to the boundary segments It is important to mention that ‘‘ks(x)’’ is treated as an RSF in the simulation. In consequence the dependent variables ‘‘h(x)’’ and ‘‘q(x)’’ also become RSFs. Now in Eq. (8) flux is a function of the conductivity and the partial of the hydraulic head; which can be obtained as a function of the pressure ‘‘p’’, density ‘‘q’’, the gravitational acceleration ‘‘g’’, viscosity ‘‘l’’, etc., as is shown in Eqs. (10) and (11). Although, some problems are simplified if it is consider a single phase and uncompressible fluid. These models are called black oil models. pðxÞ hðxÞ ¼ ð10Þ þ x3 ; qg kðxÞqg : ð11Þ k s ðxÞ ¼ l The fluid distribution is allowed along ‘‘x’’- and ‘‘y’’-axes but along the ‘‘z’’ axis is neglected due to the 2D model used. The direction of the fluid flow and the application of these equations are shown in Fig. 6. Here a defined volume is shown where new fluid is injected into a single one face and distributed to the neighbor cells as is indicated; nevertheless inside the core sample this condition changes according to the number of neighbors and their properties. The solution of the equations is subject to the following assumptions. The fluid used is incompressible; this simplify the Eqs. (4)–(9) eliminating the influence of the density in the fluid flow (q = constant). Gravity forces are neglected. Capillary does affect diffusivity or conductivity because the connecting areas between the nodes are the same ðAI;J ¼ A neighbor ¼ dz dy ¼ dz dxÞ. There is no chemical reaction between the core sample and the fluid. Profile of fluid distribution is a function of the sample properties (/, k and S). Injected fluid is the same than the resident with constant properties. Fluid injection begins at (t = 0). The easiest case for distribution is shown in Fig. 6; the nearest neighbors are in the following positions: ðI 1; J Þ; ðI þ 1; J Þ; ðI; J 1Þ and ðI; J þ 1Þ. Nevertheless, there are four more elements surrounding a pivoting node not only those in perpendicular directions as shown in Fig. 7. Due to this the maximum number of available neighbors for distribution change to eight. Although this number can be different if the analyzed node is placed in the border (wall) or in a corner of the sample. Finally, this number can be reduced if there are some neighbors not available
Fig. 6. Volume control for the simulation of fluid flow.
1760
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
Fig. 7. Neighbors of a pivoting element.
for transport process. This assumption is for those elements where the values of their properties are zero (/I;J ¼ 0Þ or ðk I;J ¼ 0Þ. They are considering as isolated elements or solid volumes, respectively, and they do not participate on the fluid flow distribution.
Fig. 8. Flow charts for the calculation of the fluid flow inside PM.
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1761
Although some authors [9,15] have developed models including different sizes of the analyzed cells and classified them in classes too, here it is considered that, this procedure is not absolutely necessary. Other authors [12,14] used a 1D model and these models are easy to solve but the fluid displacement is allowed in only one direction. A better approach is obtained using 2D models. The mathematical process include three different subroutines for the calculation of fluid flow in order to solve the partial Eq. (12) and the continuity equations for an uncompressible fluid where the variation of fluid volume is calculated for all the mesh at every iteration as a function of the simulating time (t + Dt) and the corresponding distribution velocities. In other words a new quantity of fluid is injected into the sample at every step time (Dt) until the receiving node is absolutely saturated.
Fig. 9. Evolution of the fluid injection inside PM.
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1762
o2 hð0Þ ðxÞ ohf ðxÞi ohð0Þ ðxÞ gðxÞ þ ¼ : ox2i oxi oxi K G ðxÞ
ð12Þ
The easiest case involves only the first algorithm shown in Fig. 8a; this process is used for filling of the nodes defined as inlets; it adds a new quantity of fluid as a function of the injection flux rate and at every step time. Here the second term of Eqs. (4)–(7) is not present; and the result is only the mass of new fluid added to the non-saturated inlet until this volume is exceeded. When the pore volume is totally saturated (Vpore = Vsat I,J) the algorithm of Fig. 8b is used so that the excess of fluid is distributed. The excess of fluid is distributed towards the non-saturated neighbors; if there is no neighbor available the subroutine research proceeds until it finds one or more connected and available nodes to the analyzed one, that is, when there is one or more neighbors available to the distribution (nnd = 1) or (nnd > 1), respectively. The search parameter will be equal to one (a = 1). When there is no one available neighbor the parameter will be increased to be more than one (a > 1) for searching the nearest, new available and connected neighbors. In the algorithm is considered that the distribution of fluid is affected by the distance between the analyzed and distributed nodes; so that the factor (fd) must be calculated using Eq. (13). This value is calculated for the number of neighbors for distribution (nnd) and is the sum of the partial permeability; at the end of algorithm (7b) the corresponding percent of this sum is used to calculate the mass of the excess fluid driven to the corresponding neighbor. 1 nnd1 fdnnd þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : I;J ¼ fdx;y 2 ½I x þ ½J y2
ð13Þ
Finally, there could be nodes defined as outlets in the mesh. The excess of fluid that arrives to them is counted and eliminated the calculation, the flow charts of the Fig. 8c shows this process. The calculation process will show the simulation until the user stops the process. This process is used to evaluate the displacement of the resident fluid. If there is no outlet in, the sample the simulation will conclude when all the pore volumes are completely saturated. Here the superscripts (t and t 1) mean the status of a variable in the process and in the last iteration. 5. Simulation The sample dimensions are (D = 2*, r = 25.4 mm) and (H = 38.1 mm), the total volume for simulation is (307.257 mm3) and the total porous volume is (12.60 mm3) the saturated condition was neglected. It was considered that the new fluid is injected by three inlets placed along the supply face using a constant supply (0.210 mm3/min). Fig. 9 shows the evolution of the fluid injection, it can be observed that the new fluid does not follow a constant path, the dispersion phenomena is present and the advance front will search the easiest way to continue filling the available neighbors. In a homogeneous medium the advance front of an injected fluid is well defined; nevertheless, the profile of the advance front cannot be easy to predict due to the chaotic characteristics inside a porous media. 6. Conclusions and comments The evolution of the fluid injection inside a PM was simulated successfully and the model can be scaled to meso or macro-scales easily or can be used as a platform to build a 3D model. Fluid displacement does not present a defined advance front, due to the heterogeneous distribution of its properties. Distribution of the fluid is a function of the neighbors participating in the transport and their properties; nevertheless, it is important to remember that fingering could be due to the conditions of mixing between two different fluids or when the involved phases are miscible with one another. The dispersion observed is smoothed due to that new fluid injected by three inlets along the radial face of the core sample; this distribution applies to the water jets over a larger PM area. As is shown at the beginning many branches are formed but after a few minutes they tend to converge and try to consolidate a single region of saturated fluid but with a dispersed front. Simulation of fluid flow in a PM using stochastic methods gives an advantage to obtain result about displacement performance without using more expensive techniques like tomography for virtual reconstructions [12] of PM. Acknowledgements The authors thank Instituto Polite´cnico Nacional (IPN-COFAA) and Instituto Mexicano del Petro´leo (IMP).
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1763
References [1] Pires Adolfo P. A splitting technique for analytical modelling of two-phase multi-component flow in porous media. J Petroleum Sci Eng 2006;51:54–67. [2] Hunt AG. Percolative transport in fractal porous media. Chaos, Solitons & Fractals 2004;19(2):309–25. [3] Hunt AG. Comment on fractal approach to hydraulic properties in unsaturated porous media. Chaos, Solitons & Fractals 2006;28:278–81. [4] Katz AJ, Thompson AH. Quantitative prediction of permeability in porous rock. Phys Rev B 1986;34:8179–81. [5] Mandelbrot B. The fractal geometry of nature. San Francisco Cal: Freeman; 1982. [6] Carman PC. Flow of gases through porous media. London: Butterworths; 1956. [7] Standnes Dag Chun. Spontaneous imbibition of water into cylindrical cores with high aspect ratio: numerical and experimental results. J Petroleum Sci Eng 2006;50:151–60. [9] Behbahani Hassan Sh et al. Simulation of counter-current imbibition in water-wet fractured reservoirs. J Petroleum Sci Eng 2006;50:21–39. [10] Andrade JS et al. Phys Rev Lett 1997;79(20):3901–4. [11] Andrade JS et al. Phys Rev Lett 1999;82(26):5249–52. [12] Datan JF. 3D reconstitution of porous media from image processing data using a multiscale percolation system. J Petroleum Sci Eng 2004;42:15–28. [13] Koch DL, Ladd AJC. Moderate Reynolds number flows through periodic and random arrays of aligned cylinders. J Fluid Mech 1997;349:31. [14] Othman MR, Numbere DT. Characterization of macro-scale heterogeneity and homogeneity of porous media employing fractal geometry. Chaos, Solitons & Fractals 2002;13(4):845–52. [15] Sahimi Muhammad. Flow and transport in porous media and fractured rock. VCH; 1995. [16] Sahimi Muhammad. Applications of percolation theory. London: Taylor & Francis; 1994. [17] Xu YF, Dong P. Fractal approach to hydraulic properties in unsaturated porous media. Chaos, Solitons & Fractals 2004;19(2):327–37. [18] Xu YF. Closure to fractal approach to hydraulic properties in unsaturated porous media. Chaos, Solitons & Fractals 2006;28:282–4. [19] Recommended practices for core analysis. 2nd ed. American Petroleum Institute (API); 1998. [20] Ramirez A, et al. Mathematical simulation of oil reservoir properties. Chaos, Solitons & Fractals 2008;38(3):778–88.
Simulation of uncompressible fluid flow through a porous media A. Ramı´rez
a,* ,
J.L. Gonza´lez a, F. Carrillo b, S. Lopez
c
a
b
Instituto Polite´cnico Nacional (SEPI-ESIQIE-IPN), Unidad Profesional Zacatenco, Laboratorio de Ana´lisis Met. (Edif. ‘‘Z’’ y Edif. ‘‘6’’ P.B.), Me´xico City, Mexico Instituto Polite´cnico Nacional (SEPI-CICATA-IPN), Unidad Altamira Tamaulipas, Me´xico, Mexico c Instituto Mexicano del Petro´leo (I.M.P.-D.F.), Me´xico, Mexico Accepted 18 June 2007
Abstract Recently, a great interest has been focused for investigations about transport phenomena in disordered systems. One of the most treated topics is fluid flow through anisotropic materials due to the importance in many industrial processes like fluid flow in filters, membranes, walls, oil reservoirs, etc. In this work is described the formulation of a 2D mathematical model to simulate the fluid flow behavior through a porous media (PM) based on the solution of the continuity equation as a function of the Darcy’s law for a percolation system; which was reproduced using computational techniques reproduced using a random distribution of the porous media properties (porosity, permeability and saturation). The model displays the filling of a partially saturated porous media with a new injected fluid showing the non-defined advance front and dispersion of fluids phenomena. Ó 2009 Published by Elsevier Ltd.
1. Introduction The analysis of the fluid flow behavior in a porous medium involves the solution of some mathematical equations that must be calculated simultaneously to solve transport conditions for an unsteady state. For engineering problems, the solution is quickly found when the problem is related for a homogeneous medium and the boundary conditions are easy to define and represent. The equation of continuity for fluids can be solved using a finite difference or finite element method on rectangular grids, whose solution is easy to program using nested loops algorithms. But the problem becomes complex when a porous medium is analyzed. Because the properties involved in the transport of fluid flow are not constant. Averages of porosity, permeability and saturation conditions are frequently taken to represent the behavior of a system of porous media especially for analysis on macro or mega-scales; nevertheless, these methods not always represent in detail the behavior of the fluid flow inside. For example, in natural oil reservoirs with heterogeneous properties, there are irregular porous regions bordered by regions with minor porosities or permeabilities, for this reason they exist in the world as natural confined volumes. They can be formed by groups of big rocks with *
Corresponding author. Tel.: +55 57 29 6000x55127. E-mail addresses: [email protected], [email protected] (A. Ramı´rez).
0960-0779/$ - see front matter Ó 2009 Published by Elsevier Ltd. doi:10.1016/j.chaos.2007.06.105
1754
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
available space between them or long granular regions that allow the residence of a fluid inside under high pressure conditions. Here, one of the biggest problems is to describe the geometry and properties of these reservoirs. The shape of a rock or little grains and terrain geometry, all of them are irregular and cannot be analyzed or described using simple Euclidean or analytical geometry. This means that they cannot be represented as spheres, squares, ellipsoid or any regular geometrical form. Furthermore, in the real world, natural reservoirs are absolutely saturated and due to their irregular geographical position there is an irregular pressure distribution; this pressure is very high and produces oil inside springs up to the surface. A very popular kind of models to simulate fluid flow in anisotropic materials [14,5,11] was developed at the beginning of the 1940s with fractals and stochastic theories and it is called the ‘‘random walker’’ or ‘‘drunken walker’’. In this fluid flow through the porous medium is simulated considering a relative probability for following a non-deterministic path. Some of these models have been used to represent some natural phenomenon like the cracking propagation, the randomly way of rays, etc. [4,5]. The first step to represent fluid flow systems is to define the available spaces of the medium, the outlets, the inlets and their injection rates, so that the fluid injected into the system is named the new fluid or the invasion fluid by some authors [15,16] and the resident fluid is called the old fluid or the defender fluid. The new fluid begins to conquer the nearest points of the mesh and push out the resident fluid. These models were the platform to develop the percolation theory and are based on the employment of Monte Carlo methods to predict the probable fluid path [2,3]. These models have been studied mainly by oil companies to predict the oil behavior and evaluate and improve the efficiency of their recovery techniques. Another way to simulate fluid flow is to create a computational representation of the porous media and then calculate the fluid displacement by solving fluid flow equations as a function of the medium properties. This process is illustrated in this work. Here, a section of the porous media can be reproduced taken as a base the methods of Monte Carlo and random numbers generators included in logical algorithms in order to build a group of regions with an interconnected porosity in an original non-homogeneous medium, and these regions represent the possible channels where the fluid flow will be driven. These models as a function of the scale used represent the tortuosity of the paths and the hydraulic conductivity as a parameter to calculate fluid mobility (k), permeability of the medium, etc [6,7]. Random walk models can also be programmed to display a group of cracks or a group of rocks with interfacial spaces. The points where the cracks are joined or separated in many branches will be taken as the fluid way. The same consideration can be taken for the points that define the rocks borders; here the equation of conservation is used to solve the distribution of the fluid as a function of the new fluid received and then the excess of fluid on saturated capacity will be redistributed [10–13].
2. Mathematical model for PM Here a cylindrical core of sand berea is taken as a sample as shown in Fig. 1 with closed faces. A squared mesh of a finite difference method is used to represent a split of the porous media. The dimensions and volume that represent every node are the same and remain constant during simulation. This process is repeated for assigning the three most important properties of the porous media to every node (porosity, permeability and saturation). The 2D mathematical model used to mesh the sample is described in detail in a previous work [20]. The algorithms include a subroutine for
Fig. 1. Volume of the cores sample representing the mathematical model.
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1755
Table 1 Probabilities for defining porosity of the sample Class
Position I=1
I = nx
J=1
J = ny
1 2 3 4 5
40 20 15 15 10
40 20 15 15 10
20 20 20 20 20
20 20 20 20 20
/effe 4.00 2.00 6.00 8.00 0.00
Table 2 Probabilities for defining permeability of the sample Class
Position
1 2
K
I=1
I = nx
J=1
J = ny
50 50
50 50
50 50
50 50
5 2
reading and calculating the distribution of properties. The information for the distribution is shown in Tables 1 and 2 and was taken from real data of sandstone berea cores characterization. A reproduction of the core sample is generated using an algorithm based on a stochastic assignation value subroutine. Here every value of porosity is classified in classes (c) as is done for populations in a statistical analysis, according to their presence frequency in the real core sample. Then a probability to be present in PM is defined for the perpendicular faces. So that, the appearing probability of the porosity (PcI,J) is calculated using Eq. (1) as a function of the element position (I,J) in the sample. The process is illustrated in Fig. 2 and it is repeated for all the defined classes; the probability (Pc) is obtained by the sum of the particular probabilities due to this (PcI,J = 100). Instead the probabilities are compared with a generated random number (Z), and if the value of (Z) is found in the range (Pc 1 < Z < Pc) the value of the class (property) is assigned to the node I ðP /I¼1 þ P/I¼nx Þ þ nyJ ðP /J ¼1 þ P /J¼nx Þ nx : ð1Þ PcI;J ¼ 2 Some authors [1,9,17,18] have assumed empty core samples in their models according to their particular cases to be solved. Here, we consider a partially saturated medium with a non-equal probability for distribution of the saturation classes along the longitudinal axis ‘‘x’’ but constant for the transversal axis ‘‘y’’ in contrast with the profiles obtained for porosity and permeability that remain constant along both axes as is shown, respectively, in Fig. 3a–c. A general saturation profile as a function of the maximum pore volume defined (Vmax pore) is shown in Fig. 4a. Eq. (2) is used to compare the saturation values in order to display a new saturation profile as a function of the maximum pore volume defined. Additionally, isolated or non-accessible regions of the sample can be represented with blocked nodes. These isolated elements can be easily identified in a new equivalent representation as is shown in Fig. 4b. Isolated elements are commonly used to watch all the interconnected pores (effective porosity), V sat tot I;J ¼
V sat1;J : V max pore
ð2Þ
This process has a mathematical base that is described next. It is well known that fluid flow though porous medium is related with medium spatial variation; in consequence its behavior (displacement and advance front) cannot be accurately predicted using only deterministic methods. This work is focused on describing fluid flow in a complex porous media that produce non-defined advanced fronts as a consequence of a heterogeneous distribution of the properties involved. Uncertainties in the values of the porous media properties can be reduced but it can never be entirely eliminated; the theory of stochastic process provides a natural method for evaluating them, and a tool to simulate these cases using computational techniques. In the present work the properties of the porous media were treated as random space functions (RSFs), whose spatial correlations can be quantified by probabilistic distributions as was mentioned [20]. A better explanation is given next.
1756
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
Fig. 2. Flow chart of the subroutine to define PM properties.
Fig. 3. Profiles and distribution curves for core sample properties.
A random variable ‘‘U’’ can be defined by a set ‘‘X’’ of possible values (possible values of the properties) and a function of probability distribution designed by ‘‘P’’, which is defined on a collection of subsets of ‘‘X’’ denoted by ‘‘S’’; which means that it must be closed under countable unions and finite intersections, due to this ‘‘P’’ must satisfy certain
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1757
Fig. 4. Profiles (a) saturation and (b) isolated elements of PM.
conditions. For example if ‘‘P(X) = 0’’ this denotes an empty set, if (S1 S2) then P(S1) 6 P(S2) this implicates that the event ‘‘U 2 S’’. The probability that ‘‘U’’ lies between two different values that limits a range ‘‘u’’ and ‘‘v’’ can be written as follows: ‘‘P(u 6 x 6 v) = F(v) 6 F(u)’’. Here ‘‘F’’ is the cumulative density function; now if the difference between ‘‘u’’ and ‘‘v’’ is infinitesimal ‘‘du’’ the probability that ‘‘U’’ lies between ‘‘u’’ and ‘‘du’’ is ‘‘F(u + du) F(u) = dF(u)’’ and its probability of density function will be P(u) = dF(u)/du. Finally, the probability that ‘‘U’’ lies between an interval ½a; b can be written in an integral form as is shown: P ða 6 U 6 bÞ ¼
Z
b
P ðuÞ du ¼ F ðbÞ F ðaÞ:
ð3Þ
a
A procedure based on the comparison of the probabilities to lie between a pair of values was used to define and simulated the porous media properties in a simple computational algorithm. The node definition has two main purposes. The first is to store the particular properties of everyone and the second is to establish the starting conditions for the simulation of fluid flow. In other words it must be defined which nodes will work as inlets, pores and outlets and what are their corresponding injection rates. It means that it must be defined what nodes will participate in the fluid flow and what will be their roles during distribution of exceed of fluid.
3. Characterization of PM It is important to mention that there are two kinds of porosities. The absolute porosity (/abs) is the total empty volume of a material and the effective porosity (/effe), which is formed by all connected pore spaces from an initial point to an ending point in a PM. Both properties were characterized in this work and the procedures are the following. For the (/abs), samples of sandstone berea were cut, cleaned and polished resulting in thin plates which were prepared by introducing a colored resin into their porous structure, previously evacuated, and the porous regions are identified as the geometrical zones where the resin color remains. Thus, only the opened porous spaces are accessible to fluid. After this the samples were re-polished; then a set of five photographs of one face of five samples were snapped using optical microscopy. Finally, the relation between colored area with resin (Aresin) and the area of the photograph (Atotal) represents the measure of the (/abs) and it was counted using a software image analyzer (see Table 3). For the (/effe), a procedure according to the Archimides immersion method established on the recommended practices for core analysis [19] shown in Fig. 5 was followed to analyze and classify the samples. In consequence of the calculations it is possible to confirm (/abs > /effe), because there is a volume of unconnected spaces that contributes to the (/abs) but not to the (/effe). The results of the analysis are shown in Table 4. The average was (/effe = 4.10). This information is necessary to define the porosity data used to generate the simulation shown in Fig. 3a and eventually to solve the fluid flow equations.
Table 3 Probabilities for defining saturation of the sample Class
1 2
Position
S
I=1
I = nx
J=1
J = ny
80 20
20 80
80 20
20 80
80 90
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1758
Fig. 5. Archimides mercury immersion apparatus.
Table 4 Presence of (/effe) and (/abs) inside PM Class
(/abs)
(%)
(/effe)
(%)
1 2 3 4 5
10 12 14 16 18
10.00 26.67 33.33 20.00 10.00
0 2 4 6 8
10.00 20.00 40.00 15.00 15.00
4. Mathematical model for fluid flow displacement Here we used the finite difference method to solve Eqs. (4)–(7) for every node used in the discretization process at every iteration. The left terms of these equations represent the fluid flow and movement of the new fluid in a new node and the right terms the fluid flow and movement of the resident fluid ejected from the node (when their pore volume is absolutely saturated). Eqs. (4) and (5) are the continuity equations that can be solved as a function of the quantity of movement Eqs. (6) and (7) o qu dy dz ¼ qu þ ðquÞ dx dy dz; ð4Þ ox o ð5Þ qv dx dz ¼ qv þ ðqvÞ dy dx dz; oy o quV dy dz ¼ quV þ ðquV Þ dx dy dz; ð6Þ ox o qvV dx dz ¼ qvV þ ðqvV Þ dy dx dz: ð7Þ oy Continuity equations can be written as Eq. (8) and solved as a function of Darcy’s law (9). Here ‘‘k’’ is the conductivity that can be established as a function of the permeability and subjected to the conditions shown next. r qðxÞ ¼ gðxÞ; qi ðxÞ ¼ k s
ohðxÞ ; oxi
hðxÞ ¼ H B ðxÞ ðaÞ x 2 CD ; qðxÞ nðxÞ ¼ QðxÞ ðbÞ x 2 CN :
ð8Þ ð9Þ
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1759
These equations can be solved using a deterministic set for dependent variables ‘‘h(x)’’ and ‘‘q(x)’’ if ‘‘k’’ is known at every point of a grid, where the notations used means: x is the vector of space coordinates q is the flux discharged at point (x) $ is the grad operator with respect to (x) h(x) is the the hydraulic head kS(x) is the saturated hydraulic conductivity n(x) = (n1, . . ., nd)T is an outward unit vector normal to the boundary segments It is important to mention that ‘‘ks(x)’’ is treated as an RSF in the simulation. In consequence the dependent variables ‘‘h(x)’’ and ‘‘q(x)’’ also become RSFs. Now in Eq. (8) flux is a function of the conductivity and the partial of the hydraulic head; which can be obtained as a function of the pressure ‘‘p’’, density ‘‘q’’, the gravitational acceleration ‘‘g’’, viscosity ‘‘l’’, etc., as is shown in Eqs. (10) and (11). Although, some problems are simplified if it is consider a single phase and uncompressible fluid. These models are called black oil models. pðxÞ hðxÞ ¼ ð10Þ þ x3 ; qg kðxÞqg : ð11Þ k s ðxÞ ¼ l The fluid distribution is allowed along ‘‘x’’- and ‘‘y’’-axes but along the ‘‘z’’ axis is neglected due to the 2D model used. The direction of the fluid flow and the application of these equations are shown in Fig. 6. Here a defined volume is shown where new fluid is injected into a single one face and distributed to the neighbor cells as is indicated; nevertheless inside the core sample this condition changes according to the number of neighbors and their properties. The solution of the equations is subject to the following assumptions. The fluid used is incompressible; this simplify the Eqs. (4)–(9) eliminating the influence of the density in the fluid flow (q = constant). Gravity forces are neglected. Capillary does affect diffusivity or conductivity because the connecting areas between the nodes are the same ðAI;J ¼ A neighbor ¼ dz dy ¼ dz dxÞ. There is no chemical reaction between the core sample and the fluid. Profile of fluid distribution is a function of the sample properties (/, k and S). Injected fluid is the same than the resident with constant properties. Fluid injection begins at (t = 0). The easiest case for distribution is shown in Fig. 6; the nearest neighbors are in the following positions: ðI 1; J Þ; ðI þ 1; J Þ; ðI; J 1Þ and ðI; J þ 1Þ. Nevertheless, there are four more elements surrounding a pivoting node not only those in perpendicular directions as shown in Fig. 7. Due to this the maximum number of available neighbors for distribution change to eight. Although this number can be different if the analyzed node is placed in the border (wall) or in a corner of the sample. Finally, this number can be reduced if there are some neighbors not available
Fig. 6. Volume control for the simulation of fluid flow.
1760
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
Fig. 7. Neighbors of a pivoting element.
for transport process. This assumption is for those elements where the values of their properties are zero (/I;J ¼ 0Þ or ðk I;J ¼ 0Þ. They are considering as isolated elements or solid volumes, respectively, and they do not participate on the fluid flow distribution.
Fig. 8. Flow charts for the calculation of the fluid flow inside PM.
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1761
Although some authors [9,15] have developed models including different sizes of the analyzed cells and classified them in classes too, here it is considered that, this procedure is not absolutely necessary. Other authors [12,14] used a 1D model and these models are easy to solve but the fluid displacement is allowed in only one direction. A better approach is obtained using 2D models. The mathematical process include three different subroutines for the calculation of fluid flow in order to solve the partial Eq. (12) and the continuity equations for an uncompressible fluid where the variation of fluid volume is calculated for all the mesh at every iteration as a function of the simulating time (t + Dt) and the corresponding distribution velocities. In other words a new quantity of fluid is injected into the sample at every step time (Dt) until the receiving node is absolutely saturated.
Fig. 9. Evolution of the fluid injection inside PM.
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1762
o2 hð0Þ ðxÞ ohf ðxÞi ohð0Þ ðxÞ gðxÞ þ ¼ : ox2i oxi oxi K G ðxÞ
ð12Þ
The easiest case involves only the first algorithm shown in Fig. 8a; this process is used for filling of the nodes defined as inlets; it adds a new quantity of fluid as a function of the injection flux rate and at every step time. Here the second term of Eqs. (4)–(7) is not present; and the result is only the mass of new fluid added to the non-saturated inlet until this volume is exceeded. When the pore volume is totally saturated (Vpore = Vsat I,J) the algorithm of Fig. 8b is used so that the excess of fluid is distributed. The excess of fluid is distributed towards the non-saturated neighbors; if there is no neighbor available the subroutine research proceeds until it finds one or more connected and available nodes to the analyzed one, that is, when there is one or more neighbors available to the distribution (nnd = 1) or (nnd > 1), respectively. The search parameter will be equal to one (a = 1). When there is no one available neighbor the parameter will be increased to be more than one (a > 1) for searching the nearest, new available and connected neighbors. In the algorithm is considered that the distribution of fluid is affected by the distance between the analyzed and distributed nodes; so that the factor (fd) must be calculated using Eq. (13). This value is calculated for the number of neighbors for distribution (nnd) and is the sum of the partial permeability; at the end of algorithm (7b) the corresponding percent of this sum is used to calculate the mass of the excess fluid driven to the corresponding neighbor. 1 nnd1 fdnnd þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : I;J ¼ fdx;y 2 ½I x þ ½J y2
ð13Þ
Finally, there could be nodes defined as outlets in the mesh. The excess of fluid that arrives to them is counted and eliminated the calculation, the flow charts of the Fig. 8c shows this process. The calculation process will show the simulation until the user stops the process. This process is used to evaluate the displacement of the resident fluid. If there is no outlet in, the sample the simulation will conclude when all the pore volumes are completely saturated. Here the superscripts (t and t 1) mean the status of a variable in the process and in the last iteration. 5. Simulation The sample dimensions are (D = 2*, r = 25.4 mm) and (H = 38.1 mm), the total volume for simulation is (307.257 mm3) and the total porous volume is (12.60 mm3) the saturated condition was neglected. It was considered that the new fluid is injected by three inlets placed along the supply face using a constant supply (0.210 mm3/min). Fig. 9 shows the evolution of the fluid injection, it can be observed that the new fluid does not follow a constant path, the dispersion phenomena is present and the advance front will search the easiest way to continue filling the available neighbors. In a homogeneous medium the advance front of an injected fluid is well defined; nevertheless, the profile of the advance front cannot be easy to predict due to the chaotic characteristics inside a porous media. 6. Conclusions and comments The evolution of the fluid injection inside a PM was simulated successfully and the model can be scaled to meso or macro-scales easily or can be used as a platform to build a 3D model. Fluid displacement does not present a defined advance front, due to the heterogeneous distribution of its properties. Distribution of the fluid is a function of the neighbors participating in the transport and their properties; nevertheless, it is important to remember that fingering could be due to the conditions of mixing between two different fluids or when the involved phases are miscible with one another. The dispersion observed is smoothed due to that new fluid injected by three inlets along the radial face of the core sample; this distribution applies to the water jets over a larger PM area. As is shown at the beginning many branches are formed but after a few minutes they tend to converge and try to consolidate a single region of saturated fluid but with a dispersed front. Simulation of fluid flow in a PM using stochastic methods gives an advantage to obtain result about displacement performance without using more expensive techniques like tomography for virtual reconstructions [12] of PM. Acknowledgements The authors thank Instituto Polite´cnico Nacional (IPN-COFAA) and Instituto Mexicano del Petro´leo (IMP).
A. Ramı´rez et al. / Chaos, Solitons and Fractals 39 (2009) 1753–1763
1763
References [1] Pires Adolfo P. A splitting technique for analytical modelling of two-phase multi-component flow in porous media. J Petroleum Sci Eng 2006;51:54–67. [2] Hunt AG. Percolative transport in fractal porous media. Chaos, Solitons & Fractals 2004;19(2):309–25. [3] Hunt AG. Comment on fractal approach to hydraulic properties in unsaturated porous media. Chaos, Solitons & Fractals 2006;28:278–81. [4] Katz AJ, Thompson AH. Quantitative prediction of permeability in porous rock. Phys Rev B 1986;34:8179–81. [5] Mandelbrot B. The fractal geometry of nature. San Francisco Cal: Freeman; 1982. [6] Carman PC. Flow of gases through porous media. London: Butterworths; 1956. [7] Standnes Dag Chun. Spontaneous imbibition of water into cylindrical cores with high aspect ratio: numerical and experimental results. J Petroleum Sci Eng 2006;50:151–60. [9] Behbahani Hassan Sh et al. Simulation of counter-current imbibition in water-wet fractured reservoirs. J Petroleum Sci Eng 2006;50:21–39. [10] Andrade JS et al. Phys Rev Lett 1997;79(20):3901–4. [11] Andrade JS et al. Phys Rev Lett 1999;82(26):5249–52. [12] Datan JF. 3D reconstitution of porous media from image processing data using a multiscale percolation system. J Petroleum Sci Eng 2004;42:15–28. [13] Koch DL, Ladd AJC. Moderate Reynolds number flows through periodic and random arrays of aligned cylinders. J Fluid Mech 1997;349:31. [14] Othman MR, Numbere DT. Characterization of macro-scale heterogeneity and homogeneity of porous media employing fractal geometry. Chaos, Solitons & Fractals 2002;13(4):845–52. [15] Sahimi Muhammad. Flow and transport in porous media and fractured rock. VCH; 1995. [16] Sahimi Muhammad. Applications of percolation theory. London: Taylor & Francis; 1994. [17] Xu YF, Dong P. Fractal approach to hydraulic properties in unsaturated porous media. Chaos, Solitons & Fractals 2004;19(2):327–37. [18] Xu YF. Closure to fractal approach to hydraulic properties in unsaturated porous media. Chaos, Solitons & Fractals 2006;28:282–4. [19] Recommended practices for core analysis. 2nd ed. American Petroleum Institute (API); 1998. [20] Ramirez A, et al. Mathematical simulation of oil reservoir properties. Chaos, Solitons & Fractals 2008;38(3):778–88.