Monte Carlo simulation of flow of fluids through porous media
Monte Carlo simulation of flow of fluids through porous media Semant Jain, Madhav Acharya1, Sandeep Gupta, Ashok N. Bhaskarwar * Department of Chemical Engineering, Indian Institute of Technology, Hauz Khas, New Delhi 110016, India Received 30 July 1999; received in revised form 26 August 2002; accepted 17 September 2002
Abstract
This simulation employs Monte Carlo technique for studying fluid flow through a porous medium in the capillary regime. The medium has been modelled as a 2 or 3-dimensional network of elements, some of which are randomly closed to the fluid flow. Dijkstra's algorithm has been employed to identify the least-resistance pathway, which is instrumental in determining the minimum pressure required to achieve break-through across the network. At higher pressures, network resistance has been calculated by determining the manner in which the cluster forms and by accounting for the nature of flowpaths. The simulation yields a linear relationship between the pressure applied across the network and flowrate showing similarity to Darcy's law. Polynominal fitting of the data on the fraction of openable pores open as dependent on pressure applied across the network has been carried out and the coefficients determined. Keywords: Monte Carlo; Simulation; Porous media; Fluid flow; Darcy's law; Scaling
1. Introduction
Monte Carlo technique is a method of computer simulation of a system with many degrees of freedom. It makes use of random numbers to numerically generate probability distributions, which might otherwise not be explicitly known since the considered systems are so complex (Binder, 1979). Monte Carlo simulation provides a good comparison between data from experiments on real systems to those from the model. It is used in areas like simulation of thermodynamic properties of fluids, crystal growth, combustion of coal particles etc. Its classical application includes evaluation of multiple integrals in statistical mechanics.
The term 'porous media' encompasses a wide variety of contacting devices such as packed towers, sand beds and substances like limestone rock, filter paper and catalytic particles. It is desirable to classify the porous media according to the types of pore spaces they contain. A proposed classification was by dividing the pore spaces into voids, capillaries and force spaces (Manegold, 1937). Void spaces are characterized by the fact that walls have little or no effect on hydrodynamic properties in the interior; in capillaries, the walls do affect the hydrodynamics but do not bring the molecular structure of the fluid into evidence; and in force spaces, the molecular structure of the fluid is of considerable importance. This work concentrates on pores of the size of capillaries. Flow through a porous medium requires a description of both the medium and the flow. A porous medium can be represented as an extremely complicated network of channels, including those containing obstructions and dead ends too (Bernsdorf, Brenner & Durst, 2000). The distribution of channels is obtained from assumed statistical descriptions. The pores in the network can be interconnected or non-connected, depending on whether they are a part of a continuous network of pores that exists within the medium or not. Put another
386
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
Nomenclature
Small letters
A n r
h Capital letters Pa Pc Pcr -* max
Q Qi Qmax
Ri
fraction of pores that have been blocked number of flowpaths present capillary radius size of the grid along one dimension or side applied pressure (on inlet face of network) capillary pressure break-through pressure maximum pressure total flowrate flowrate of ith flowpath maximum flowrate through the network radius of curvature of meniscus resistance of ith flowpath
Greek symbols a
surface tension contact angle
way, any channel classified as 'interconnected' will ultimately be filled by fluid flowing through the medium, while 'non-interconnected' elements will remain devoid of any fluid flow. The fraction of interconnected channels gives an indication of the accessible porosity of the medium. To describe the flow through a porous medium, we also need to specify two parameters—applied pressure across the network and the flowrate (i.e. the net amount of fluid passing through the network per unit time). We have divided the paper in six main sections. The theory section has an overview of percolation theory, capillary pressure and pore structure models. The section following that describes the approach we used to develop the simulation code. It covers the algorithm used to identify the least resistant pathway, how the status of a pore changes from 'closed' initially to 'open' to finally 'part of a flowpath', the manner in which subsequent flowpaths are identified and finally the equations used to compute the flowrate through the medium. After briefly mentioning our constraints, we describe our results and finally present the conclusions.
2. Theoretical background
2.1. Percolation theory There are many physical phenomena in which a fluid spreads randomly through a medium, e.g. it may be a solute diffusing through a solvent, molecules penetrating a porous solid, or electrons migrating over an atomic lattice. Besides the random mechanism, external forces
may also govern the process, as in the case of water flowing through limestone under gravity. According to the nature of the problem, the random mechanism can be attributed to the fluid or the medium. The former falls under the category of diffusion processes. A typical example is the motion of one molecule in a gas as it undergoes collisions with other molecules. In case of a dilute gas, each collision event is totally random as it is not influenced by other collisions that have occurred in the past. In other words, the medium has no 'memory' of its past history. Also, the medium (which is essentially the molecules), is continuously varying after each collision and so is not invariant in time. The other, relatively less common, is known as percolation. In percolation processes (such as a fluid soaking into a porous medium), there is a distinction between the fluid particles and the scattering medium. This medium, although it varies in random fashion from point to point, is invariant in time. Thus 'memory' effects cannot be neglected as in diffusion, and the random scattering of the particles of the fluid must be treated as being an inherent property of the medium. This difference between percolation and diffusion can be mathematically understood through the 1-dimensional Polya walk. The medium is described as a set of points placed at equal intervals along a straight line and the particles of fluid can move in steps of unit length in either direction with equal probability. In the case of diffusion, the points constitute the fluid as well as the medium and so they can move in random fashion without any constraints. In the corresponding percolation process, the points of the medium are assigned a direction to the left or right with equal probability. A particle entering the medium moves in accordance with
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
387
b
o Fig. 1. Particles 'a' and 'b' are trapped in the medium due to orientation of arrows.
O • O
o o
Interconnected elements Oiled with fluid Non interconnected elements
O • O
Interconnected elements not filled by fluid Fig. 3. Intermediate stage of percolation in 2-dimensional network.
Fig. 2. Percolation pathway found in the network. Pores part of pathway have been shown connected with a dashed line. Black filled circles are inter-connected pores. Open circles are blocked pores. Arrows indicate direction of motion of a particle at a point in the medium.
the arrows at each point and so the medium plays the active role. As can be seen from Fig. 1, a particle can get trapped and be forced to oscillate indefinitely between two points (whereas this does not occur in the case of diffusion). In such a case, no percolation path can be struck between two ends of a line. This simple 1-dimensional illustration can be extended further to 2 and 3-dimensional networks of elements, where some elements may also be blocked off to fluid flow. Consider a 2-dimensional matrix through which a path has to be found (Fig. 2). If all elements are open to fluid flow, then they form part of a single 'percolation path' for the entire medium. As elements are blocked off, the size of the percolation path is reduced, but it still connects both ends of the network. At a certain critical fraction of closed elements, the percolation path ceases to exist. It has been observed that this fraction is 0.41 for 2-dimensional networks and 0.69 for 3-dimensional ones (Efros, 1982) (Fig. 3).
2.2. Theory of capillary pressure Consider the hydrostatics of two immiscible fluids or phases that exist simultaneously in a porous medium (Greenkorn, 1983; De Weist, 1969). In general, one phase will wet the solid. The entrance of one fluid into a small pore against the other fluid is opposed by surface tension forces between the wetting fluid and the pore walls (Scheidegger, 1963; Muskat, 1982). The result is that a certain pressure differential in the displacing phase versus the displaced phase will have to be produced to maintain equilibrium. This pressure is called the capillary pressure. In a single capillary, the curvature Rc of the interface gives rise to the pressure differential equal to Pc = 2s=Rc
(1)
The radius of curvature of the meniscus is equal to Rc = r=cos f
(2)
So that for a single circular capillary Pc = 2s cos f=r
(3)
As can be seen from the above expression, the surface tension force is inversely related to capillary radius. Hence, capillary pressure can be regarded as the resistance offered by a capillary to the flow of fluid through it—the larger the capillary radius, the lower the
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
resistance. In this simulation, the resistance of an element is dimensionless and so surface tension values are irrelevant, i.e. the analysis is system non-specific. 2.3. Pore structure models This simulation uses a simplified model of porous medium. The model consists of a network of elements that represent cylindrical capillary tubes of different diameters and equal lengths. A single pore, then, is a series of elements placed one after another and so incorporates the effect of varying diameter along its length. Our model has all pores of same length. The effect of 'arrows' (as in Polya walk) is obtained by the random generation and assignment of resistances (i.e. through assigning radii values) to elements of the medium. The concept of least resistance is used to determine the percolation path of the fluid. If pressure is applied to a fluid-filled porous medium, or to the fluid at the entrance to a capillary system, the fluid will penetrate those capillaries whose capillary pressure is lower than the applied pressure. In other words, the largest-diameter capillaries would be filled first and at increasing pressures, the smaller capillaries would get filled. This is referred to as the concept of 'least resistance'. In actual porous media, the pores can be fully filled, partially filled or be completely empty. Although the simulation assumes a pore under consideration as being either fully filled or completely empty, it is possible to model a partially filled pore as a combination of two adjacent pores—one being fully-filled and the other completely empty. In this simulation, the pressure applied across the network is incremented by a very small amount at every iteration. This approach is validated by experimental observations where the effect of hysteresis was diminished or eliminated by carrying out the experiment sufficiently slowly (Dullien, 1992). Varying pore geometry in a porous medium essentially implies a pore of a greater or smaller 'capillary pressure' than a corresponding pore of uniform diameter. This implies that this pore requires a greater or smaller pressure for break-through. Such effects have been accounted for by selecting the minimum and maximum pore radii at the beginning of the simulation. The random selection of a radius value between the above limits for a pore incorporates all media effects that can affect the resistance of a pore (Dullien). Thus, though this simplified model does not account for the presence of films, 'permanent hysteresis' (Dullien, 1992) or different connectivities of the pore elements with their neighbours (Cordero, Rojas & Riccardo, 2001) in the system, it incorporates most features of porous media that would have a strong bearing on the flowrate.
3. Simulation 3.1. Initialization The simulation requires a random generation of resistances of the elements that mimic the porous medium. The parameters of relevance are the size of the network or lattice and the desired phase fraction of blocked elements (Monteagudo, Rajagopal & Lage, 2002). The resistances are generated over a range of values as would exist in a porous medium, but the actual distributions as reported in literature have not been used. At a given instant, each pore is either 'closed', 'open' or 'part of a flowpath'. To start with, all nodes are marked 'closed'. It should be noted that there might be pore(s) having finite resistance but are surrounded completely by 'blocked-off pores. As these pores cannot be reached from the entry face, their status would remain unchanged from 'closed' throughout the duration of the simulation. As the pressure is increased and the fluid begins to percolate in the medium, pores that are filled are marked 'open'. For a pore's status to become 'part of a flowpath', it must become a part of either a dependent or an independent flowpath. It is possible for pores to remain 'open' and yet not 'part of a flowpath' because these pores could be dead ends or a sequence of pores that have as yet not succeeded in forming a flowpath. 3.2. Least resistant pathway The simulation employs Dijkstra's algorithm to determine the least resistant pathway. Since both pressure and resistance are dimensionless, the resistance value of the least resistant pathway can be equated to give the 'break-through' or the minimum pressure required to cause the first flowpath to appear. The simulation has an entry face for the fluid to enter into the network which can be regarded as a 'singlesource' for the algorithm. By definition, the resistances of all pores are positive implying positive edge weights for the graph. The concept of edge weight is equivalent to the resistance of a pore and is stored in the pore itself. The simulation maintains a priority queue that contains all the pores whose source distance is yet to be finalized. All pores present in the queue are sorted in ascending order of their source distance. To begin with, all pores are assigned a source distance value as infinite. Then the source distance of the pores on the entry face is equated to the resistance of the pore. This is followed by repeatedly selecting a pore with the least source distance in the priority queue and relaxing all pores that are connected to it. The process of relaxation of the source distance of a pore involves updating the current value of the source distance with the sum of the resistance of the
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
pore and the source distance of the neighbouring pore (that has just been popped from the queue), if the value of the latter is smaller (Cormen, Leiserson & Rivest, 2000). After all the pores have been assigned the minimum possible source distance value, the source distance of all pores on the exit face is compared and the one having the least value is selected. This is the value of resistance of the 'least resistant pathway' or the 'break-through' pressure required across the network. For a pathway to exist between the entry and exit face of the network, the value of 'break-through' pressure must be finite. If a pathway has been found, the elements are stored in order of their appearance in the percolation path and consequently, the order in which they will be filled up by the fluid. 3.3. Opening of a node The next step is to carry out a depth-first search across all pores having finite resistance which have opened up. A pore is deemed to have 'opened-up' when its capillary pressure is smaller than the effective pressure available at that pore. Effective pressure at a pore is the difference of the cumulative resistance of all pores that precede the current pore from the applied pressure across the network. The recursion for a pore ends when the capillary resistance is greater than the effective pressure or all adjacent pores having finite resistance have been explored. The pressure is incremented in steps until all the pores that have finite resistance and are reachable have 'opened-up'. Although the pressure is incremented slowly, at a given increment there can be multiple pores that 'open-up'. This is similar to the morphological approach used to study fractal dimensions (Hilpert & Miller, 2001). 3.4. Additional flowpath determination The process of repeatedly increasing applied pressure would eventually lead to other flowpaths opening up. A new flowpath comes into existence when an independent or dependent flowpath comes into existence. An 'independent flowpath' consists of a series of pores that have 'opened-up' from the entry to the exit face, while a 'dependent flowpath' is a sequence of pores from the entry face to a pore already part of an existing flowpath. The resistance offered by an 'independent flowpath' is the sum of the resistances of all pores that form the flowpath from the entry to the exit face. For a 'dependent flowpath', the simplest case would involve just one sequence of pores from the entry face to a pore already part of an existing flowpath. Here, the resistance of the newly opened flowpath equals the cumulative resistance of all pores that have just recently
389
'opened-up' and the 'flow resistance' of the pore already a part of some existing flowpath. 'Flow resistance' of a pore is defined as the cumulative resistance of all the pores starting from the pore under consideration to a pore on the exit face. Other than the above mentioned case, there is a distinct possibility of the existence of multiple series of pores that have opened up from a pore inside the network and joined existing flowpaths. In order to take into account the branching of flow from such pores, the resistance of each branch is computed separately by linear addition of the resistances of the pores that have 'opened-up' (starting at the pore from which the branching begins till a pore which is an immediate neighbour of a pore whose status is 'part of a flowpath') with the 'flow-resistance' of the pore which is a 'part of a flowpath'. Thus, the 'flow resistance' of the pore from which the branching begins is the parallel summation of the resistances of all the branches from that pore onwards. The effective resistance of the 'dependent flowpath' is thus the linear addition of the resistances of all pores from the entry face to the pore from which branching begins and the 'flow resistance' of that pore. 3.5. Flowrate
In the case of an actual porous medium, the physical quantity 'flowrate' is defined only when the fluid actually exits from the pores at the end opposite to the one at which it entered (assuming 1-dimensional percolation). The simulated network, in reality, has several exit points, each of which has its own flowrate. The flowrate through the network is defined as the cumulative flowrate from all the pores on the exit face. In this simulation, mass and volume conservation have been assumed to hold. The pores do not rupture in the pressure range being studied. Since fluid enters from one face and leaves from the opposite face, the cumulative flowrate could be alternatively defined as the sum of the flowrates entering the network through the pores on the entry face. When the first flowpath is obtained (using Dijkstra's algorithm), the net flowrate through it (and in this case, through the network) is zero. As the inlet pressure is increased, a new sequence of pores starts getting filled by the fluid, which may result in another flowpath joining one of the existing flowpath or flowpaths emerging from the exit face. All flowpaths have the same inlet and outlet pressures at any given time and vary only in their individual resistances. Thus, the net pressure driving force across a particular path, rather than the inlet pressure is taken for calculation of the flowrate. The flowrate through a path is then computed by dividing the driving force by the path resistance. Pi = P,~Ri
(4)
390
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
150x150 Blockage 20 % 11 10 09 a. O o 0,7
S 0-5
a.
° 04
y = 4E-20X6 + 5E-16x5 - 3E-12x4 • 7E-09x3 - 9Efl6x2 * 0.0064x - 0.6577 R2 = 0.978
£ 02 01 00 500
1000
1500
2000
2500
3000
3500
Pressure
Fig. 4. Fraction of openable pores open vs. applied pressure for a 150 x 150 grid at 20 % blockage. Continuous line represents the least-squares averaged trend line.
150x150 Blockage 20% 90000 80000 y = 21.098x- 353.29 70000 -
R2 = 0357
60000 S
50000 -
j2
40000 30000 20000 10000 0 500
1500
1000
2000
2500
3000
3500
Pressure
Fig. 5. Flowrate vs. applied pressure for a 150 x 150 grid at 20% blockage. Continuous line represents the least squares averaged trend line.
'~
P
'
'
The total flowrate is then the sum of individual flowrates and can be expressed as
(6) The flowrate is calculated each time a pressure increment is made and also a new pore 'opens up'. The pressure after which no more pores open up is called 'maximum pressure'. Once 'maximum pressure' is
attained, the matrix is regenerated and the process repeated a few dozen times over. 3.6. Constraints The random nature of the trials causes variation in the values of 'break-through' and 'maximum' pressures. Both these pressures vary with matrix size and fraction of pores blocked-off. The choice of matrix size and blocked fraction is dictated by two factors—computer memory size and critical phase fraction values. The memory factor limits
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
391
80x80 Blockage 30%
c
0.9 -
Ope
1.0 •
0.8 -
w
n
o 0.7 a.
mw
0.6 -
I
0.5 •
Q.
w
o
racti
c 0.4 o 0.3 -
U.
0.2 -
y = ^E-Wx 6 * 1E-14x5 4E-11X4 * 5EJ)8x3 3Efl5x2 • 0.0126x 0.7981 R2 = 0.9206
0.1 -
-
0.0 -
200
400
600
800
1000
1200
1400
1600
1800
Pressure Fig. 6. Fraction of openable pores open vs. Applied pressure for 80 x 80 grid at 30% blockage. Continuous line represents the least-square averaged trend line.
80x80 Blockage 30% 18000 16000 14000 12000 "S
10000 -
a: O
8000 -
y=4.7191x -20.149 R2 = 0.7591
LL.
6000 4000 2000 0 0
400
600
800
1000
1200
1400
1600
1800
Pressure Fig. 7. Flowrate vs. applied pressure for 80 x 80 grid at 30 % blockage. Continuous line represents the least-squares averaged trend line.
the maximum network size that can be simulated but there is also a certain minimum value below which the simulation is unable to generate statistically significant data to mimic the flow properly. At these 'finite' sizes, the randomness associated with Monte Carlo simulation cannot be used effectively and errors result.
4. Results
In this work, we have a pore throat diameter distribution ('pore size distribution') which is generated by a pseudorandom number generator according to a uniform distribution between the minimum and max-
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
392
70x70 Blockage 40% 0.9 1.8 a. 0.7 -
oM a>
0.6 -
j>
0.5 -
enab
o
04 -
n
l-ra ion
o
0.3 -
<j
0.2 y = -6E-19x 6 + 4E-15x 5 -1E-11 x 4 + 1 E-08x3 -1 E-05x2 + 0.0048x - 0.1933 0.1 -
R2 = 0.5226
0 -
200
400
600
800
1000
1200
1400
1600
Pressure
Fig. 8. Fraction of openable pores open vs. applied pressure for 70 x 70 grid at 40% blockage. Continuous line represents the least-squares averaged trend line.
70x70 Blockage 40% 8000
7000 6000 -
y = 0.9645x-83.534 [^ = 02412
5000 4000 3000 2000 1000 0 200
400
600
800
1000
1200
1400
1600
Pressure
Fig. 9. Flowrate vs. applied pressure for a 70 x 70 grid at 40% blockage. Continuous line represents the least-squares averaged trend line.
imum values. These values were assumed to be 1 and 10, respectively. The pseudo-random numbers are used for assigning values to the required throat diameters, so that those are uncorrelated, i.e. the size of one throat is independent of the size of any other throat. Porous media may have polymodal and/or spatially correlated pore size distributions, that can also be handled by the model, but these are not considered in this work. The flowrate vs. applied pressure relationships obtained from the model for different sets of parameters show a linear dependence similar to Darcy's law (Figs. 5,
7, 9 and 11). In these plots, the values of the y -intercept must be negative. It is indicative of the fact that only at a certain finite positive pressure does the flow through the porous medium begin to take place. The values of parameters that have been estimated using linear leastsquares technique, show a consistency in the predicted values regardless of the matrix size for a specified blockage fraction (Tables 1-5). Due to simulation time constraints we could not carry out runs at grid sizes greater than 200 for 2-dimensional networks. However looking at the graphs of the dominant para-
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
393
18x18x18 Blockage 4 0 % 1.1
pen
1.0 0.9 -
o
0.8 -
in a>
0.7 -
pen ble
o a.
0.6 -
o ro
0.4 -
fl10
0.3 -
u.
0.2 -
10
0.5 -
y= -6E-16X8 + 1E-12x5 - 1EJ)9x4 + 7E-07x3 -0.0002x2 + 0.0243x - 0.1639 R2 = 0.9738
0.1 -
nn 200
100
300
400
500
600
700
Pressure Fig. 10. Fraction of openable pores open vs. applied pressure for 18 x 18 x 18 grid at 40% blockage. Continuous line represents the least-squares averaged trend line.
18x18x18 Blockage40% 1 iUUU
10000-
8000-
£ 6000-
y=14.005x 114.95 R2 = 0.9703
4000-
2000-
0-
a
100
200
300
400
500
B00
700
Pressure Fig. 11. Flowrate vs. applied pressure for a 18 x 18 x 18 grid at 40% blockage. Continuous line represents the least-squares averaged trend line.
meter vs. grid size (Fig. 12) and cumulative flowrate vs. grid size at a specific applied pressure (Fig. 13), we believe that an infinite medium can be simulated satisfactorily in the 200-300 grid size range for 2dimensional networks. The flowrate is dependent on the number of flowpaths and that in turn on the accessible porosity of the medium. In our view, fraction of openable pores open (FOPO) is a good measure of the accessible porosity. To calculate FOPO, we noted the actual number of pores that have 'opened-up' (APO) while incrementing the pressure and then divided it by the number of pores that
are not 'blocked-off’. For 2-dimensional systems; Openable Pores = sg2(1 —fh) For 3-dimensional systems; Openable Pores = ^ ( 1 —fb) FOPO = APO=Openable Pores
(7)
(8) (9)
On plotting FOPO vs. applied pressure for different configurations, i.e. grid sizes and blockage fractions, we note an initial linear rising trend, attaining an asymptotic value subsequently (Figs. 4, 6, 8 and 10). This trend is
Table 1 Values of coefficients for fraction of openable pores open and flowrate for 2 and 3-dimensional grids at blockage of 20% Dimension
Size
Fraction of openable pores open
Flowrate R2
Coefficients
3
P5
P4
P3
P2
P1
P0
40 50 60 70 80 90 100 120 150 175 190 200
-3.00E-16 -8.00E-17 -2.00E-17 -9.00E-18 -4.00E-18 -2.00E-18 -8.00E-19 -2.00E-19 -4.00E-20 -5.00E-21 1.00E-20 1.00E-20
1.00E-12 3.OOE-13 1.00E-13 5.00E-14 2.00E-14 1.00E-14 6.00E-15 2.00E-15 5.00E-16 1.00E-16 -2.00E-16 -1.00E-16
-1.00E-09 -4.00E-10 -2.00E-10 -1.00E-10 -6.00E-11 -3.00E-11 -2.00E-11 -8.00E-12 -3.00E-12 -8.00E-13 6.00E-13 5.00E-13
7.00E-07 3.00E-07 2.00E-07 1.00E-07 7.00E-08 5.00E-08 3.00E-08 2.00E-08 7.00E-09 3.00E-09 -6.00E-10 -6.00E-10
-2.00E-04 -1.00E-04 -9.00E-05 -6.00E-05 -5.00E-05 -3.00E-05 -3.00E-05 -2.00E-05 -9.00E-06 -5.00E-06 -1.00E-06 -8.00E-07
3.28E-02 2.52E-02 2.08E-02 1.72E-02 1.46E-02 1.25E-02 1.11E-02 8.60E-02 6.40E-02 4.60E-02 2.50E-02 2.20E-02
-9.49E-02 -8.47E-01 -8.66E-01 -8.27E-01 -8.03E-01 -7.73E-01 -7.65E-01 -7.08E-01 -6.58E-01 -5.15E-01 -2.16E-01 -1.85E-01
10 11 12 13 14 15 16 17 18 19 20
-2.00E-13 -3.00E-14 -2.00E-14 -1.00E-15 -5.00E-15 -2.00E-15 -1.00E-15 -4.00E-16 -3.00E-17 1.00E-16 8.00E-17
2.00E-10 5.00E-11 3.00E-11 2.00E-11 8.00E-12 5.00E-12 3.00E-12 1.00E-12 3.OOE-13 1.00E-13 -8.00E-14
-8.00E-08 -2.00E-08 -2.00E-08 -1.00E-08 -6.00E-09 -4.00E-09 -2.00E-09 -1.00E-09 -5.00E-10 -9.00E-11 -6.00E-11
2.00E-05 6.00E-06 4.00E-06 3.00E-06 2.00E-06 1.00E-06 1.00E-06 6.00E-07 4.00E-07 2.00E-07 1.00E-07
-1.50E-03 -8.00E-04 -6.00E-04 -5.00E-04 -4.00E-04 -3.00E-04 -2.00E-04 -2.00E-04 -1.00E-04 -7.00E-05 -6.00E-05
7.16E-02 4.88E-02 4.28E-02 3.75E-02 3.14E-02 2.89E-02 2.46E-02 2.10E-02 1.77E-02 1.45E-02 1.31E-02
-2.30E-01 -1.03E-01 -8.75E-02 -8.07E-02 -2.83E-02 -4.22E-02 -2.17E-02 1.70E-03 3.88E-02 -4.50E-03 8.50E-03
Q
R2 ain el
2
P6
Coefficients P1
P0
0.9089 0.9306 0.9441 0.9539 0.9606 0.9641 0.9679 0.9746 0.9780 0.9777 0.9727 0.9748
4.89 5.98 7.81 9.28 10.51 11.65 12.88 16.79 21.10 23.58 26.98 28.20
/29.30 /40.76 /47.56 /41.79 /48.02 /81.98 /77.17 /44.28 /353.29 /259.32 /590.17 /120.83
4863.60 5935.94 7760.15 9235.21 10456.98 11567.02 12797.83 16744.72 20744.71 23324.68 26384.83 28082.17
0.8195 0.7727 0.8829 0.9014 0.9080 0.8840 0.9263 0.9195 0.9570 0.9112 0.9728 0.9568
0.9821 0.9496 0.9644 0.9680 0.9735 0.9752 0.9799 0.9805 0.9831 0.9637 0.9670
11.45 13.97 17.06 20.25 23.28 27.24 30.44 35.76 39.38 44.23 50.05
/99.17 /77.61 /96.06 /110.19 /85.03 /143.68 /164.09 /187.90 /212.43 /252.98 /283.51
11349.83 13892.40 16959.95 20136.81 23297.97 27097.32 30276.91 35573.10 39162.57 43973.02 49761.49
0.9570 0.9681 0.9729 0.9679 0.9680 0.9765 0.9786 0.9842 0.9805 0.9725 0.9684
R2 represents the regression coefficient values. Q represents the value of flowrate at an applied pressure of 1000 units. FOPO (Pa) = P6i 3 a+P5i 3 a+P4Pa + P3Pa P1Pa + P0, Q = P1 x 1000+P0.
P0, Flowrate (P a ) =
6.
!
e
a. 2
1 eg
1 ^ j
o ^~
Table 2 Values of coefficients for fraction of openable pores open and flowrate for 2- and 3-dimensional grids at blockage of 25 % Dimension
Size
Fraction of openable pores open
Flowrate
Coefficients
Coefficients
P6
P5
P4
P2
P1
P0
P1
P0
Q
R2
40 50 60 70 80 90 100 120 150 175 200
-3.00E-16 -7.00E-17 -2.00E-17 -8.00E-18 -4.00E-18 -2.00E-18 -8.00E-19 -2.00E-19 -4.00E-20 +9.00E-21 9.00E-21
8.00E-13 3.OOE-13 1.00E-13 4.00E-14 2.00E-14 1.00E-14 6.00E-15 2.00E-15 5.00E-16 -6.00E-17 -1.00E-16
-1.00E-09 -4.00E-10 -2.00E-10 -9.00E-11 -5.00E-11 -3.00E-11 -2.00E-11 -8.00E-12 -3.00E-12 -4.00E-14 4.00E-13
6.00E-07 3.00E-07 2.00E-07 1.00E-07 7.00E-08 4.00E-08 3.00E-08 2.00E-08 7.00E-09 1.00E-09 -2.00E-10
-2.00E-04 -1.00E-04 -8.00E-05 -6.00E-05 -4.00E-05 -3.00E-05 -3.00E-05 -2.00E-05 -9.00E-06 -3.00E-06 -1.00E-06
3.05E-02 2.46E-02 2.01E-02 1.70E-02 1.49E-02 1.28E-02 1.11E-02 8.60E-03 6.40E-03 3.8OE-O3 2.50E-03
-8.67E-01 -9.09E-01 -8.97E-01 -8.84E-01 + 8.95E-01 + 8.39E-01 -8.45E-01 -7.74E-01 -7.25E-01 -4.23E-01 -2.67E-01
0.8789 0.9162 0.9325 0.9444 0.9540 0.9618 0.9626 0.9710 0.9761 0.9708 0.9715
3.37 4.60 5.33 6.22 7.58 8.93 9.06 12.04 14.72 16.81 21.83
/23.16 /23.86 /34.30 /49.19 /38.44 /123.30 /54.21 /141.63 /67.14 /1.86 /103.21
3344.14 4575.55 5298.41 6174.31 7540.87 8811.10 9007.69 11902.57 14647.86 16810.14 21725.79
0.7737 0.7250 0.8331 0.7869 0.8264 0.8281 0.8748 0.9463 0.8777 0.9140 0.9471
10 11 12 13 14 15 16 17 18 19 20
-2.00E-13 -4.00E-14 -2.00E-14 -1.00E-14 -5.00E-15 -3.00E-15 -1.00E-15 -6.00E-16 -2.00E-16 1.00E-16 2.00E-16
2.00E-10 5.00E-11 3.00E-11 2.00E-11 9.00E-12 5.00E-12 3.00E-12 1.00E-12 7.00E-13 -6.00E-14 -2.00E-13
-7.00E-08 -2.00E-08 -2.00E-08 -1.00E-08 -6.00E-09 -4.00E-09 -2.00E-09 -1.00E-09 -8.00E-10 -2.00E-10 2.00E-11
1.00E-05 6.00E-06 4.00E-06 3.00E-06 2.00E-06 1.00E-06 1.00E-06 1.00E-07 5.00E-07 2.00E-07 1.00E-07
-1.50E-03 -8.00E-04 -6.00E-04 -5.00E-04 -4.00E-04 -3.00E-04 -2.00E-04 -2.00E-04 -1.00E-04 -8.00E-05 -6.00E-05
7.05E-02 5.07E-02 4.34E-02 3.83E-02 3.43E-02 2.91E-02 2.60E-02 2.23E-02 1.97E-02 1.56E-02 1.38E-02
-2.24E-01 -1.47E-01 -1.12E-01 -1.06E-01 -1.03E-01 -6.88E-02 -5.82E-02 -2.82E-02 -1.02E-02 -3.82E-02 -1.75E-02
0.9696 0.9550 0.9602 0.9673 0.9694 0.9726 0.9774 0.9813 0.9816 0.9596 0.9642
10.30 11.34 13.95 16.49 19.91 21.35 25.06 27.20 30.77 34.90 38.75
/46.54 /61.61 /79.27 /94.49 /111.98 /127.20 /146.22 /160.49 /183.64 /216.07 /244.84
10253.47 11282.40 13873.73 16392.51 19796.02 21222.80 24916.78 27038.51 30585.36 34680.93 38504.16
0.9596 0.9627 0.9470 0.9704 0.9684 0.9694 0.9800 0.9834 0.9806 0.9673 0.9660
R2 represents the regression coefficient values. Q represents the value of flowrate at an applied pressure of 1000 units. POPO (Pa) = P1Pa + P0, Q = Pl x 1000+P0.
, Flowrate (Pa) =
I3. 0
Table 3 Values of coefficients for fraction of openable pores open and flowrate for 2- and 3-dimensional grids at blockage of 30 % Size
Fraction of openable pores open
Flowrate R2
Coefficients
3
P5
P4
P3
P2
P1
P0
40 50 60 70 80 90 100 120 150 175 190 200
-2.00E-16 -5.00E-17 -9.00E-18 -5.00E-18 -2.00E-18 -9.00E-19 -4.00E-19 -1.00E-19 -2.00E-20 -1.00E-21 1.00E-20 1.00E-20
6.00E-13 2.00E-13 5.00E-14 3.00E-14 1.00E-14 7.00E-15 3.00E-15 1.00E-15 2.00E-16 6.00E-17 -1.00E-16 -1.00E-16
-7.00E-10 -3.00E-10 -9.00E-11 -6.00E-11 -4.00E-11 -2.00E-11 -1.00E-11 -6.00E-12 -2.00E-12 -6.00E-13 5.00E-13 6.00E-13
5.00E-07 2.00E-07 1.00E-07 7.00E-08 5.00E-08 3.00E-08 2.00E-08 1.00E-08 5.00E-09 2.00E-09 -5.00E-10 -9.00E-10
-2.00E-04 -1.00E-04 -5.00E-05 -4.00E-05 -3.00E-05 -3.00E-05 -2.00E-05 -1.00E-05 -8.00E-06 -5.00E-06 -1.00E-06 -2.00E-07
2.64E-02 2.14E-02 1.54E-02 1.44E-02 1.26E-02 1.09E-02 9.60E-03 8.00E-03 6.00E-03 4.40E-03 2.40E-03 1.80E-03
-7.83E-01 -8.51E-01 -7.23E-01 -8.11E-01 -7.98E-01 -8.05E-01 -7.90E-01 -7.87E-01 -7.42E-01 -6.10E-01 -2.68E-02 -1.83E-01
10 11 12 13 14 15 16 17 18 19 20
-5.00E-14 -3.00E-14 -2.00E-14 -1.00E-14 -5.00E-15 -3.00E-15 -2.00E-15 -3.00E-12 -3.00E-16 8.00E-18 1.00E-16
6.00E-11 4.00E-11 3.00E-11 2.00E-11 9.00E-12 5.00E-12 4.00E-12 6.00E-09 9.00E-13 2.00E-13 -6.00E-14
-3.00E-08 -2.00E-08 -2.00E-08 -1.00E-08 -6.00E-09 -4.00E-09 -3.00E-09 -6.00E-06 -1.00E-09 -4.00E-10 -1.00E-10
7.00E-06 6.00E-06 4.00E-06 3.00E-06 2.00E-06 2.00E-06 1.00E-06 +2.80E-03 5.00E-07 3.00E-07 2.00E-07
-8.00E-04 -8.00E-04 -6.00E-04 -5.00E-04 -4.00E-04 -3.00E-04 -3.OOE-O3 -6.85E-01 -2.00E-04 -1.00E-04 -8.00E-05
4.90E-02 5.16E-02 4.44E-02 4.01E-02 3.43E-02 3.06E-02 2.80E-02 8.28E+01 2.11E-02 1.77E-02 1.55E-02
-7.78E-02 -1.93E-01 -1.51E-01 -1.52E-01 -1.22E-01 -1.12E-01 -1.13E-01 -2.70E+02 -5.63E-01 -9.66E-02 -6.76E-02
Q
R2 ain el
2
P6
Coefficients P1
P0
0.8487 0.8784 0.8944 0.9227 0.9203 0.9383 0.9534 0.9572 0.9679 0.9728 0.9657 0.9643
2.39 3.12 3.43 4.04 4.72 5.43 6.90 7.84 10.44 13.02 12.42 13.69
/22.77 /19.96 /28.82 /34.01 /20.15 /40.26 /48.15 /39.59 /126.99 /34.94 /71.08 /144.91
2363.33 3104.14 3396.98 4004.39 4698.95 5392.14 6855.05 7804.61 10316.01 12980.06 12352.92 13544.09
0.6176 0.6818 0.7054 0.8041 0.7591 0.8179 0.8026 0.8617 0.8970 0.9100 0.8820 0.8988
0.9693 0.9671 0.9558 0.9599 0.9885 0.9732 0.9723 0.9777 0.9628 0.9581 0.9613
7.63 9.69 10.75 12.98 14.71 16.80 18.70 21.56 24.47 27.45 30.88
/68.05 /78.84 /72.81 /86.76 /98.58 /114.86 /123.03 /141.11 /137.24 /191.75 /213.58
7558.55 9615.77 10680.19 12896.24 14614.42 16684.14 18576.97 21415.89 24334.76 27262.25 30666.42
0.9357 0.9626 0.9530 0.9738 0.9621 0.9601 0.9706 0.9740 0.9765 0.9839 0.9707
R2 represents the regression coefficient values. Q represents the value of flowrate at an applied pressure of 1000 units. FOPO (PB) = P6P6a + />5i3a + P4i 3a + P3i 3 a +P2i 3 a+Pli 3 a + P0, Flowrate (Pa) = PI/* / P0, 0 = P1 x 1000+ P0
6. |
5a a.
3
1
1 •00
Dimension
^~
Table 4 Values of coefficients for Fraction of Openable Pores Open and Flowrate for 2- and 3-dimensional grids at blockage of 35 % Dimension
Size
Fraction of openable pores open
Flowrate
Coefficients
2
40 50 60 70 80 90 100 120 150 175 200
3
10 11 12 13 14 15 16 17 18 19 20
R
P6
P5
P4
P3
P2
P1
P0
-2.00E-16 -9.00E-18 -4.00E-17 -7.00E-19 -1.00E-18 2.00E-19 -4.00E-19 1.00E-19 6.00E-20 2.00E-20 1.00E-20
5.00E-13 5.00E-14 5.00E-14 7.00E-15 9.00E-15 -5.00E-16 3.00E-15 -7.00E-16 -5.00E-16 -3.00E-16 -2.00E-16
-6.00E-10 -9.00E-11 1.00E-10 -2.00E-11 -2.00E-11 -1.00E-12 -9.00E-12 1.00E-12 2.00E-12 1.00E-12 8.00E-13
3.00E-07 9.00E-08 1.00E-07 3.00E-08 3.00E-08 6.00E-09 1.00E-08 3.00E-10 -2.00E-09 -2.00E-09 -2.00E-09
-1.00E-04 -5.00E-05 -5.00E-05 -2.00E-05 -3.00E-05 -9.00E-06 -1.00E-05 -4.00E-06 -2.00E-07 1.00E-06 1.00E-06
1.95E-02 1.31E-02 1.97E-02 9.70E-03 1.00E-02 5.60E-03 6.70E-03 4.00E-03 2.20E-03 9.00E-04 9.00E-04
-5.43E-01 -4.75E-01 7.44E-02 -5.72E-01 -6.96E-01 -3.58E-01 -5.38E-01 -3.82E-01 -2.43E-01 -5.20E-03 -6.94E-02
-1.00E-13 -2.00E-14 -2.00E-14 -1.00E-14 -6.00E-15 -3.00E-15 -2.00E-15 -1.00E-15 -6.00E-16 -2.00E-16 8.00E-17
1.00E-10 3.00E-11 3.00E-11 2.00E-11 1.00E-11 6.00E-12 4.00E-12 2.00E-12 1.00E-12 6.00E-13 4.00E-15
-6.00E-08 -2.00E-08 -1.00E-08 -1.00E-08 -7.00E-09 -4.00E-09 -3.00E-09 -2.00E-09 -1.00E-09 -7.00E-10 -2.00E-10
1.00E-05 5.00E-06 4.00E-06 3.00E-06 2.00E-06 2.00E-06 1.00E-06 9.00E-07 7.00E-07 4.00E-07 2.00E-07
-1.3OE-O3 -7.00E-04 -6.00E-04 -5.00E-04 -4.00E-04 -3.00E-04 -3.00E-04 -2.00E-04 -2.00E-04 -1.00E-04 -9.00E-05
6.77E-02 4.52E-02 4.47E-02 4. HE-02 3.63E-02 3.22E-02 2.90E-02 2.57E-02 2.43E-02 2.03E-02 1.68E-02
-2.74E-01 -1.61E-01 -1.83E-01 -1.94E-01 -1.78E-01 -1.59E-01 -1.47E-01 -1.22E-01 -1.64E+01 -1.71E-01 -1.11E-01
2
Coefficients
Q
R2
P1
P0
0.6144 0.7336 0.7235 0.8182 0.8293 0.8645 0.8276 0.8983 0.9486 0.8808 0.9425
1.30 1.80 2.05 2.32 3.11 4.22 4.38 4.97 5.95 7.86 /425.30
/21.16 /16.92 /41.64 /20.71 /22.08 /19.90 /108.20 /129.83 /83.84 /483.51 8163.10
1275.94 1784.38 2010.56 2300.49 3088.02 4196.30 4268.90 4842.57 5861.36 7372.69 8163.1
0.4851 0.5058 0.5494 0.5519 0.7216 0.7665 0.6773 0.7234 0.7084 0.8720 0.8574
0.9703 0.9570 0.9499 0.9511 0.9645 0.9680 0.9690 0.9729 0.9755 0.9558 0.9565
5.62 6.63 8.18 9.54 10.97 12.78 15.49 16.74 14.01 20.71 22.93
/49.53 /41.61 /60.93 /72.61 /82.84 /92.88 /116.15 /134.24 /114.95 /162.10 /182.16
5572.47 6591.19 8120.87 9470.10 10887.16 12691.12 15377.85 16608.76 13890.05 20543.90 22747.84
0.9240 0.9507 0.9574 0.9595 0.9451 0.9546 0.9662 0.9692 0.9836 0.9772 0.9680
a.
R2 represents the regression coefficient values. Q represents the value of flowrate at an applied pressure of 1000 units. POPO (Pa) = P6Pa + P5JPa+P4JPa+P3JPa + IP2Pa 2 VlP, ,+P0, Flowrate (Pa) = PlPo / P0, 0 = P1 x 1000+ P0.
s. o o
12 ! i
3 2. eg
1
Table 5 Values of coefficients for Fraction of Openable Pores Open and Flowrate for 2- and 3-dimensional grids at blockage of 40 % Dimension
Size
Fraction of openable pores open
Flowrate
2
3
R
Coefficients P1
R2
Q
P6
P5
P4
P3
P2
P1
P0
P0
40 50 60 70 80 90 100 120 150 175
-4.00E-17 -4.00E-18 -2.00E-17 -6.00E-19 -4.00E-18 8.00E-19 -7.00E-19 3.00E-19 2.00E-19 1.00E-20
1.00E-13 6.00E-15 7.00E-14 4.00E-15 -2.00E-14 -5.00E-15 5.00E-15 -2.00E-15 -2.00E-15 -1.00E-16
-1.00E-10 9.00E-13 -1.00E-10 -1.00E-11 4.00E-11 1.00E-11 -1.00E-11 5.00E-12 6.00E-12 7.00E-13
1.00E-07 -5.00E-09 1.00E-07 1.00E-08 -4.00E-08 -1.00E-08 2.00E-08 -5.00E-09 -1.00E-08 -2.00E-09
-4.00E-05 5.00E-07 6.00E-05 1.00E-05 2.00E-05 5.00E-06 -1.00E-05 3.00E-07 1.00E-05 3.00E-06
8.40E-03 1.50E-03 1.48E-02 4.80E-03 -2.00E-03 1.00E-04 4.80E-03 2.40E-03 -6.60E-03 -7.20E-03
-7.93E-02 3.15E-01 -8.78E-01 -1.93E-01 3.63E-01 1.93E-01 -3.47E-01 -3.62E-01 1.59E+00 5.66E-01
0.4145 0.3538 0.4960 0.5226 0.5690 0.6605 0.5954 0.6712 0.6274 0.6346
0.53 1.11 0.93 0.96 1.35 2.42 1.78 3.07 3.85 5.01
/4.57 /48.56 /39.81 /83.53 /90.05 /97.95 /65.86 /871.17 /153.03 /983.12
525.53 1060.84 894.19 880.97 1262.75 2317.45 1709.44 2197.03 3697.67 4029.18
0.1420 0.4811 0.3329 0.2412 0.3091 0.3922 0.3082 0.3891 0.5664 0.5165
10 11 12 13 14 15 16 17 18 19 20
-5.00E-14 -3.00E-14 -2.00E-14 -1.00E-14 -6.00E-15 -3.00E-15 -2.00E-15 -1.00E-15 -6.00E-16 -3.00E-16 -7.00E-17
6.00E-11 5.00E-11 3.00E-11 2.00E-11 1.00E-11 6.00E-12 4.00E-12 3.00E-12 1.00E-12 8.00E-13 3.OOE-13
-3.00E-08 -2.00E-08 -2.00E-08 -1.00E-08 -7.00E-09 -5.00E-09 -3.00E-09 -2.00E-09 -1.00E-09 -9.00E-10 -5.00E-10
7.00E-06 6.00E-06 4.00E-06 3.00E-06 2.00E-06 2.00E-06 1.00E-06 1.00E-06 7.00E-07 5.00E-07 3.00E-07
-9.00E-04 -8.00E-04 -7.00E-04 -5.00E-04 -4.00E-04 -4.00E-04 -3.OOE-O3 -2.00E-04 -2.00E-04 -1.00E-04 -1.00E-04
5.27E-02 5.28E-02 4.64E-02 4.19E-02 3.78E-02 3.40E-02 3.OOE-O3 2.76E-02 2.43E-02 2.17E-02 1.91E-02
-1.81E-01 -2.68E-01 -2.34E-01 -2.45E-01 -2.35E-01 -2.27E-01 -1.96E-01 -1.94E-01 -1.64E-01 -2.19E-01 -1.85E-01
0.9378 0.9569 0.9445 0.9559 0.9605 0.9638 0.9664 0.9715 0.9738 0.9470 0.9520
3.90 5.20 5.83 7.19 8.34 9.60 11.03 12.25 14.01 16.06 17.28
/24.87 /46.21 /49.80 /59.08 /70.91 /80.99 /91.04 /101.36 /114.95 /136.46 /150.29
3870.13 5154.49 5782.70 7134.32 8266.79 9519.11 10934.97 12148.64 13890.05 15924.54 17129.71
0.9471 0.9255 0.9363 0.9373 0.9439 0.9582 0.9712 0.9659 0.9703 0.9639 0.9755
R2 represents the regression coefficient values. Q represents the value of flowrate at an applied pressure of 1000 units. POPO (Pa) = P6i3a + P5i 3 a +P4,f 1 a >P3,P a + P2 P1P a + P0, Q = Pl x 1000 + P0.
Jain
Coefficients
2
6.
, Flowrate (Pa) =
!
e
a. 2
11 . eg
1 ^ j
o ^~
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
399
Flow Rate (Slope) 100.0
4
• • •
10.0 -
•
•
•
K
•
•
n
A
A
x
•
o to
A
•
•
A
A
x
x
x
0
0
O
0
O
0
X
y.
n
x
•
• s
t
0
X
A
A X
A
A
n
A
O
0
•
o
*
1.0 -
•
•
n 20%
A
25% * 30% o 35% 140%
0.1 0
25
50
75
100
125
150
175
200
225
Grid Size Fig. 12. Slope of flowrate vs. grid size at different blockage fractions. Simulations with grid size 40 and above rate for 2-dimensional grids, while stimulations smaller than this size were for 3-dimensional grids.
FlowRate(Pa= 1000) JUUIJU
s
•
w • CgQ
10000 -
•
*
•
•
A
•A •••:
D
D
A
• A
*
•
.
n
A
y.
A
0
0
O y.
y
0
O
O
•
0
x
o
o
•
A
•
A
Q
X
1000 -
•A
• •
*
*
100 -
10 -
° 20% A 25% x 30% °35% • 40% 125
50
75
100
125
150
175
200
225
Grid Size Fig. 13. Flowrate vs. grid size at different blockage fractions at an applied pressure of 1000 units. Simulations with grid size 40 and above were for 2dimensional grids while simulations smaller than this size were for 3-dimensional grids.
expected and holds irrespective of the grid size or the fraction of pores blocked. This is similar to the observed trend when calculated cumulative mercury intrusion volume was plotted against applied pressure (Bryntesson, 2002). When the least-squares analysis is conducted and seven parameters of the polynomial estimated, it is observed that the parameter values gradually attain asymptotic values on increasing the grid size (Tables 1 -
5). This observation is in line with the fact that an ideal simulation would have all grid dimensions as infinite. At low blocked fractions, there are few chances of pores with finite resistance being surrounded by 'blocked-off pores. Thus all pores having finite resistance can be considered 'openable'. The asymptotic value of FOPO attained in all such cases is near 1 (Figs. 4 and 6). As the blockage fraction increases, the
400
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
probability of finite- resistance pores being surrounded by 'blocked-off pores is higher and thus, a greater variation is seen in the asymptotic values of FOPO. This observation is most stark at blockage values within 5% of the critical value (Fig. 8). This observation can be explained by the fact that near critical values there is a sharp decrease in the probability of finding a flowpath, the implication being that there is a greater variation in the probability of finding a flowpath for a given configuration. For this to be true, it would be natural to expect a greater variation in the asymptotic value of the openable pores. These observations have been validated by the regression analysis while computing the parameter values. The R2 trend values indicate that there is greater deviation from 1 at blockage fractions near critical values, particularly at low grid sizes (Tables 1-5). As grid size increases the R2 values improve indicating that the simulation is closer to that of an infinite medium.
5. Conclusions
This simulation is a simplified but innovative way of dealing with flow through porous media. Using the Monte Carlo technique and concepts of percolation theory, simulations have been carried out for 2 and 3dimensional models of porous media in capillary flow regime. The flowrate vs. applied pressure results have been found to be in agreement with Darcy's law. Parameter estimation of fraction of openable pores that have 'opened up' vs. applied pressure has also been shown to be consistent. More accurate results could be obtained by using actual resistance distributions of porous media. It might then be possible to predict flowrate vs. applied pressure behaviour for systems with different physical properties (e.g. surface tension and contact angle). Actual experimental data for flow in the capillary regime have not been reported in literature, where the major focus has been on statistical description of porous media. Also, experiments that reproduce the exact conditions of the simulation would be difficult to perform. Therefore, the results mentioned herein should be treated as being of a qualitative nature. The parameter values have been computed after considering several thousands of points for each configuration (matrix size and blocked fraction). Runs were concluded once the curve predicting FOPO vs. applied pressure exhibited an asymptotic behaviour. It has been our observation that until this asymptotic behaviour is observed the runs should not be concluded, otherwise the parameter values would be misleading. This is also in accordance with the random nature of the simulation
because theoretically an infinite number of runs are required before presenting any trend. Further work in this field would include extension of the simulation to the region of Hagen-Poiseulle flow, where the resistance is inversely related to the fourth power of the radius of the element. The resistance beyond the laminar flow region decreases slightly and then reaches a constant value at very high flowrates in the turbulent zone. This can be modelled as two separate flow regions, with elements having different (but constant) resistances in each region.
Acknowledgements
We are grateful to Dr M.N. Gupta, Senior Manager, Computer Services Centre, IIT, Delhi, for the generous allocation of computational time. We would also like to thank all staff members of CSC, IIT, Delhi, for their cooperation throughout this project. In particular, we would like to mention Mr. Gopal Krishen, the system programmer, and Mr. Gulshan Naveriya, the hardware engineer, for their invaluable assistance.
References Bernsdorf, J., Brenner, G., & Durst, F. (2000). Numerical analysis of the pressure drop in porous media flow with lattice Boltzmann (BGK) automata. Computer Physics Communications 129, 247255. Binder, K. (Ed.), Monte Carlo methods in statistical physics (1979). Springer-Verlag. Bryntesson, L. M. (2002). Pore network modelling of the behaviour of a solute in chromatography media: transient and steady-state diffusion properties. Journal of Chromatography A 945, 103-115. Cordero, S., Rojas, F., & Riccardo, J. L. (2001). Simulation of threedimensional porous networks. Colloids and Surfaces A: Physicochemical and Engineering Aspects 187-188, 425-438. Cormen, T. H., Leiserson, C. E., & Rivest, R. L. (2000). Introduction to algorithms. EEE: Prentice Hall of India. De Weist, R. J. M. (Ed.), Flow through porous media (1969). Academic Press. Dullien, F. A. L. (1992). Porous media—fluid transport and pore structure. Academic Press. Efros, A. L. (1982). Physics and geometry of disorder—percolation theory. Mir Publishers. Greenkorn, R. A. (1983). Flow phenomenon in porous media. Marcel Dekker. Hilpert, M., & Miller, C. T. (2001). Pore-morphology-based simulation of drainage in totally wetting porous media. Advances in Water Resources 24, 243-255. Manegold, V. E. (1937). Uber Kapillarsysteme. Kolloid-Z 80, 253. Monteagudo, J. E. P., Rajagopal, K., & Lage, P. L. C. (2002). Simulating oil flow in porous media under asphaltene deposition. Chemical Engineering Science 57, 323-337. Muskat, M. (1982). Flow of homogenous fluids through porous media.
IHRDC. Scheidegger, A. E. (1963). The physics of flow through porous media. University of Toronto Press.
Abstract
This simulation employs Monte Carlo technique for studying fluid flow through a porous medium in the capillary regime. The medium has been modelled as a 2 or 3-dimensional network of elements, some of which are randomly closed to the fluid flow. Dijkstra's algorithm has been employed to identify the least-resistance pathway, which is instrumental in determining the minimum pressure required to achieve break-through across the network. At higher pressures, network resistance has been calculated by determining the manner in which the cluster forms and by accounting for the nature of flowpaths. The simulation yields a linear relationship between the pressure applied across the network and flowrate showing similarity to Darcy's law. Polynominal fitting of the data on the fraction of openable pores open as dependent on pressure applied across the network has been carried out and the coefficients determined. Keywords: Monte Carlo; Simulation; Porous media; Fluid flow; Darcy's law; Scaling
1. Introduction
Monte Carlo technique is a method of computer simulation of a system with many degrees of freedom. It makes use of random numbers to numerically generate probability distributions, which might otherwise not be explicitly known since the considered systems are so complex (Binder, 1979). Monte Carlo simulation provides a good comparison between data from experiments on real systems to those from the model. It is used in areas like simulation of thermodynamic properties of fluids, crystal growth, combustion of coal particles etc. Its classical application includes evaluation of multiple integrals in statistical mechanics.
The term 'porous media' encompasses a wide variety of contacting devices such as packed towers, sand beds and substances like limestone rock, filter paper and catalytic particles. It is desirable to classify the porous media according to the types of pore spaces they contain. A proposed classification was by dividing the pore spaces into voids, capillaries and force spaces (Manegold, 1937). Void spaces are characterized by the fact that walls have little or no effect on hydrodynamic properties in the interior; in capillaries, the walls do affect the hydrodynamics but do not bring the molecular structure of the fluid into evidence; and in force spaces, the molecular structure of the fluid is of considerable importance. This work concentrates on pores of the size of capillaries. Flow through a porous medium requires a description of both the medium and the flow. A porous medium can be represented as an extremely complicated network of channels, including those containing obstructions and dead ends too (Bernsdorf, Brenner & Durst, 2000). The distribution of channels is obtained from assumed statistical descriptions. The pores in the network can be interconnected or non-connected, depending on whether they are a part of a continuous network of pores that exists within the medium or not. Put another
386
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
Nomenclature
Small letters
A n r
h Capital letters Pa Pc Pcr -* max
Q Qi Qmax
Ri
fraction of pores that have been blocked number of flowpaths present capillary radius size of the grid along one dimension or side applied pressure (on inlet face of network) capillary pressure break-through pressure maximum pressure total flowrate flowrate of ith flowpath maximum flowrate through the network radius of curvature of meniscus resistance of ith flowpath
Greek symbols a
surface tension contact angle
way, any channel classified as 'interconnected' will ultimately be filled by fluid flowing through the medium, while 'non-interconnected' elements will remain devoid of any fluid flow. The fraction of interconnected channels gives an indication of the accessible porosity of the medium. To describe the flow through a porous medium, we also need to specify two parameters—applied pressure across the network and the flowrate (i.e. the net amount of fluid passing through the network per unit time). We have divided the paper in six main sections. The theory section has an overview of percolation theory, capillary pressure and pore structure models. The section following that describes the approach we used to develop the simulation code. It covers the algorithm used to identify the least resistant pathway, how the status of a pore changes from 'closed' initially to 'open' to finally 'part of a flowpath', the manner in which subsequent flowpaths are identified and finally the equations used to compute the flowrate through the medium. After briefly mentioning our constraints, we describe our results and finally present the conclusions.
2. Theoretical background
2.1. Percolation theory There are many physical phenomena in which a fluid spreads randomly through a medium, e.g. it may be a solute diffusing through a solvent, molecules penetrating a porous solid, or electrons migrating over an atomic lattice. Besides the random mechanism, external forces
may also govern the process, as in the case of water flowing through limestone under gravity. According to the nature of the problem, the random mechanism can be attributed to the fluid or the medium. The former falls under the category of diffusion processes. A typical example is the motion of one molecule in a gas as it undergoes collisions with other molecules. In case of a dilute gas, each collision event is totally random as it is not influenced by other collisions that have occurred in the past. In other words, the medium has no 'memory' of its past history. Also, the medium (which is essentially the molecules), is continuously varying after each collision and so is not invariant in time. The other, relatively less common, is known as percolation. In percolation processes (such as a fluid soaking into a porous medium), there is a distinction between the fluid particles and the scattering medium. This medium, although it varies in random fashion from point to point, is invariant in time. Thus 'memory' effects cannot be neglected as in diffusion, and the random scattering of the particles of the fluid must be treated as being an inherent property of the medium. This difference between percolation and diffusion can be mathematically understood through the 1-dimensional Polya walk. The medium is described as a set of points placed at equal intervals along a straight line and the particles of fluid can move in steps of unit length in either direction with equal probability. In the case of diffusion, the points constitute the fluid as well as the medium and so they can move in random fashion without any constraints. In the corresponding percolation process, the points of the medium are assigned a direction to the left or right with equal probability. A particle entering the medium moves in accordance with
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
387
b
o Fig. 1. Particles 'a' and 'b' are trapped in the medium due to orientation of arrows.
O • O
o o
Interconnected elements Oiled with fluid Non interconnected elements
O • O
Interconnected elements not filled by fluid Fig. 3. Intermediate stage of percolation in 2-dimensional network.
Fig. 2. Percolation pathway found in the network. Pores part of pathway have been shown connected with a dashed line. Black filled circles are inter-connected pores. Open circles are blocked pores. Arrows indicate direction of motion of a particle at a point in the medium.
the arrows at each point and so the medium plays the active role. As can be seen from Fig. 1, a particle can get trapped and be forced to oscillate indefinitely between two points (whereas this does not occur in the case of diffusion). In such a case, no percolation path can be struck between two ends of a line. This simple 1-dimensional illustration can be extended further to 2 and 3-dimensional networks of elements, where some elements may also be blocked off to fluid flow. Consider a 2-dimensional matrix through which a path has to be found (Fig. 2). If all elements are open to fluid flow, then they form part of a single 'percolation path' for the entire medium. As elements are blocked off, the size of the percolation path is reduced, but it still connects both ends of the network. At a certain critical fraction of closed elements, the percolation path ceases to exist. It has been observed that this fraction is 0.41 for 2-dimensional networks and 0.69 for 3-dimensional ones (Efros, 1982) (Fig. 3).
2.2. Theory of capillary pressure Consider the hydrostatics of two immiscible fluids or phases that exist simultaneously in a porous medium (Greenkorn, 1983; De Weist, 1969). In general, one phase will wet the solid. The entrance of one fluid into a small pore against the other fluid is opposed by surface tension forces between the wetting fluid and the pore walls (Scheidegger, 1963; Muskat, 1982). The result is that a certain pressure differential in the displacing phase versus the displaced phase will have to be produced to maintain equilibrium. This pressure is called the capillary pressure. In a single capillary, the curvature Rc of the interface gives rise to the pressure differential equal to Pc = 2s=Rc
(1)
The radius of curvature of the meniscus is equal to Rc = r=cos f
(2)
So that for a single circular capillary Pc = 2s cos f=r
(3)
As can be seen from the above expression, the surface tension force is inversely related to capillary radius. Hence, capillary pressure can be regarded as the resistance offered by a capillary to the flow of fluid through it—the larger the capillary radius, the lower the
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
resistance. In this simulation, the resistance of an element is dimensionless and so surface tension values are irrelevant, i.e. the analysis is system non-specific. 2.3. Pore structure models This simulation uses a simplified model of porous medium. The model consists of a network of elements that represent cylindrical capillary tubes of different diameters and equal lengths. A single pore, then, is a series of elements placed one after another and so incorporates the effect of varying diameter along its length. Our model has all pores of same length. The effect of 'arrows' (as in Polya walk) is obtained by the random generation and assignment of resistances (i.e. through assigning radii values) to elements of the medium. The concept of least resistance is used to determine the percolation path of the fluid. If pressure is applied to a fluid-filled porous medium, or to the fluid at the entrance to a capillary system, the fluid will penetrate those capillaries whose capillary pressure is lower than the applied pressure. In other words, the largest-diameter capillaries would be filled first and at increasing pressures, the smaller capillaries would get filled. This is referred to as the concept of 'least resistance'. In actual porous media, the pores can be fully filled, partially filled or be completely empty. Although the simulation assumes a pore under consideration as being either fully filled or completely empty, it is possible to model a partially filled pore as a combination of two adjacent pores—one being fully-filled and the other completely empty. In this simulation, the pressure applied across the network is incremented by a very small amount at every iteration. This approach is validated by experimental observations where the effect of hysteresis was diminished or eliminated by carrying out the experiment sufficiently slowly (Dullien, 1992). Varying pore geometry in a porous medium essentially implies a pore of a greater or smaller 'capillary pressure' than a corresponding pore of uniform diameter. This implies that this pore requires a greater or smaller pressure for break-through. Such effects have been accounted for by selecting the minimum and maximum pore radii at the beginning of the simulation. The random selection of a radius value between the above limits for a pore incorporates all media effects that can affect the resistance of a pore (Dullien). Thus, though this simplified model does not account for the presence of films, 'permanent hysteresis' (Dullien, 1992) or different connectivities of the pore elements with their neighbours (Cordero, Rojas & Riccardo, 2001) in the system, it incorporates most features of porous media that would have a strong bearing on the flowrate.
3. Simulation 3.1. Initialization The simulation requires a random generation of resistances of the elements that mimic the porous medium. The parameters of relevance are the size of the network or lattice and the desired phase fraction of blocked elements (Monteagudo, Rajagopal & Lage, 2002). The resistances are generated over a range of values as would exist in a porous medium, but the actual distributions as reported in literature have not been used. At a given instant, each pore is either 'closed', 'open' or 'part of a flowpath'. To start with, all nodes are marked 'closed'. It should be noted that there might be pore(s) having finite resistance but are surrounded completely by 'blocked-off pores. As these pores cannot be reached from the entry face, their status would remain unchanged from 'closed' throughout the duration of the simulation. As the pressure is increased and the fluid begins to percolate in the medium, pores that are filled are marked 'open'. For a pore's status to become 'part of a flowpath', it must become a part of either a dependent or an independent flowpath. It is possible for pores to remain 'open' and yet not 'part of a flowpath' because these pores could be dead ends or a sequence of pores that have as yet not succeeded in forming a flowpath. 3.2. Least resistant pathway The simulation employs Dijkstra's algorithm to determine the least resistant pathway. Since both pressure and resistance are dimensionless, the resistance value of the least resistant pathway can be equated to give the 'break-through' or the minimum pressure required to cause the first flowpath to appear. The simulation has an entry face for the fluid to enter into the network which can be regarded as a 'singlesource' for the algorithm. By definition, the resistances of all pores are positive implying positive edge weights for the graph. The concept of edge weight is equivalent to the resistance of a pore and is stored in the pore itself. The simulation maintains a priority queue that contains all the pores whose source distance is yet to be finalized. All pores present in the queue are sorted in ascending order of their source distance. To begin with, all pores are assigned a source distance value as infinite. Then the source distance of the pores on the entry face is equated to the resistance of the pore. This is followed by repeatedly selecting a pore with the least source distance in the priority queue and relaxing all pores that are connected to it. The process of relaxation of the source distance of a pore involves updating the current value of the source distance with the sum of the resistance of the
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
pore and the source distance of the neighbouring pore (that has just been popped from the queue), if the value of the latter is smaller (Cormen, Leiserson & Rivest, 2000). After all the pores have been assigned the minimum possible source distance value, the source distance of all pores on the exit face is compared and the one having the least value is selected. This is the value of resistance of the 'least resistant pathway' or the 'break-through' pressure required across the network. For a pathway to exist between the entry and exit face of the network, the value of 'break-through' pressure must be finite. If a pathway has been found, the elements are stored in order of their appearance in the percolation path and consequently, the order in which they will be filled up by the fluid. 3.3. Opening of a node The next step is to carry out a depth-first search across all pores having finite resistance which have opened up. A pore is deemed to have 'opened-up' when its capillary pressure is smaller than the effective pressure available at that pore. Effective pressure at a pore is the difference of the cumulative resistance of all pores that precede the current pore from the applied pressure across the network. The recursion for a pore ends when the capillary resistance is greater than the effective pressure or all adjacent pores having finite resistance have been explored. The pressure is incremented in steps until all the pores that have finite resistance and are reachable have 'opened-up'. Although the pressure is incremented slowly, at a given increment there can be multiple pores that 'open-up'. This is similar to the morphological approach used to study fractal dimensions (Hilpert & Miller, 2001). 3.4. Additional flowpath determination The process of repeatedly increasing applied pressure would eventually lead to other flowpaths opening up. A new flowpath comes into existence when an independent or dependent flowpath comes into existence. An 'independent flowpath' consists of a series of pores that have 'opened-up' from the entry to the exit face, while a 'dependent flowpath' is a sequence of pores from the entry face to a pore already part of an existing flowpath. The resistance offered by an 'independent flowpath' is the sum of the resistances of all pores that form the flowpath from the entry to the exit face. For a 'dependent flowpath', the simplest case would involve just one sequence of pores from the entry face to a pore already part of an existing flowpath. Here, the resistance of the newly opened flowpath equals the cumulative resistance of all pores that have just recently
389
'opened-up' and the 'flow resistance' of the pore already a part of some existing flowpath. 'Flow resistance' of a pore is defined as the cumulative resistance of all the pores starting from the pore under consideration to a pore on the exit face. Other than the above mentioned case, there is a distinct possibility of the existence of multiple series of pores that have opened up from a pore inside the network and joined existing flowpaths. In order to take into account the branching of flow from such pores, the resistance of each branch is computed separately by linear addition of the resistances of the pores that have 'opened-up' (starting at the pore from which the branching begins till a pore which is an immediate neighbour of a pore whose status is 'part of a flowpath') with the 'flow-resistance' of the pore which is a 'part of a flowpath'. Thus, the 'flow resistance' of the pore from which the branching begins is the parallel summation of the resistances of all the branches from that pore onwards. The effective resistance of the 'dependent flowpath' is thus the linear addition of the resistances of all pores from the entry face to the pore from which branching begins and the 'flow resistance' of that pore. 3.5. Flowrate
In the case of an actual porous medium, the physical quantity 'flowrate' is defined only when the fluid actually exits from the pores at the end opposite to the one at which it entered (assuming 1-dimensional percolation). The simulated network, in reality, has several exit points, each of which has its own flowrate. The flowrate through the network is defined as the cumulative flowrate from all the pores on the exit face. In this simulation, mass and volume conservation have been assumed to hold. The pores do not rupture in the pressure range being studied. Since fluid enters from one face and leaves from the opposite face, the cumulative flowrate could be alternatively defined as the sum of the flowrates entering the network through the pores on the entry face. When the first flowpath is obtained (using Dijkstra's algorithm), the net flowrate through it (and in this case, through the network) is zero. As the inlet pressure is increased, a new sequence of pores starts getting filled by the fluid, which may result in another flowpath joining one of the existing flowpath or flowpaths emerging from the exit face. All flowpaths have the same inlet and outlet pressures at any given time and vary only in their individual resistances. Thus, the net pressure driving force across a particular path, rather than the inlet pressure is taken for calculation of the flowrate. The flowrate through a path is then computed by dividing the driving force by the path resistance. Pi = P,~Ri
(4)
390
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
150x150 Blockage 20 % 11 10 09 a. O o 0,7
S 0-5
a.
° 04
y = 4E-20X6 + 5E-16x5 - 3E-12x4 • 7E-09x3 - 9Efl6x2 * 0.0064x - 0.6577 R2 = 0.978
£ 02 01 00 500
1000
1500
2000
2500
3000
3500
Pressure
Fig. 4. Fraction of openable pores open vs. applied pressure for a 150 x 150 grid at 20 % blockage. Continuous line represents the least-squares averaged trend line.
150x150 Blockage 20% 90000 80000 y = 21.098x- 353.29 70000 -
R2 = 0357
60000 S
50000 -
j2
40000 30000 20000 10000 0 500
1500
1000
2000
2500
3000
3500
Pressure
Fig. 5. Flowrate vs. applied pressure for a 150 x 150 grid at 20% blockage. Continuous line represents the least squares averaged trend line.
'~
P
'
'
The total flowrate is then the sum of individual flowrates and can be expressed as
(6) The flowrate is calculated each time a pressure increment is made and also a new pore 'opens up'. The pressure after which no more pores open up is called 'maximum pressure'. Once 'maximum pressure' is
attained, the matrix is regenerated and the process repeated a few dozen times over. 3.6. Constraints The random nature of the trials causes variation in the values of 'break-through' and 'maximum' pressures. Both these pressures vary with matrix size and fraction of pores blocked-off. The choice of matrix size and blocked fraction is dictated by two factors—computer memory size and critical phase fraction values. The memory factor limits
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
391
80x80 Blockage 30%
c
0.9 -
Ope
1.0 •
0.8 -
w
n
o 0.7 a.
mw
0.6 -
I
0.5 •
Q.
w
o
racti
c 0.4 o 0.3 -
U.
0.2 -
y = ^E-Wx 6 * 1E-14x5 4E-11X4 * 5EJ)8x3 3Efl5x2 • 0.0126x 0.7981 R2 = 0.9206
0.1 -
-
0.0 -
200
400
600
800
1000
1200
1400
1600
1800
Pressure Fig. 6. Fraction of openable pores open vs. Applied pressure for 80 x 80 grid at 30% blockage. Continuous line represents the least-square averaged trend line.
80x80 Blockage 30% 18000 16000 14000 12000 "S
10000 -
a: O
8000 -
y=4.7191x -20.149 R2 = 0.7591
LL.
6000 4000 2000 0 0
400
600
800
1000
1200
1400
1600
1800
Pressure Fig. 7. Flowrate vs. applied pressure for 80 x 80 grid at 30 % blockage. Continuous line represents the least-squares averaged trend line.
the maximum network size that can be simulated but there is also a certain minimum value below which the simulation is unable to generate statistically significant data to mimic the flow properly. At these 'finite' sizes, the randomness associated with Monte Carlo simulation cannot be used effectively and errors result.
4. Results
In this work, we have a pore throat diameter distribution ('pore size distribution') which is generated by a pseudorandom number generator according to a uniform distribution between the minimum and max-
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
392
70x70 Blockage 40% 0.9 1.8 a. 0.7 -
oM a>
0.6 -
j>
0.5 -
enab
o
04 -
n
l-ra ion
o
0.3 -
<j
0.2 y = -6E-19x 6 + 4E-15x 5 -1E-11 x 4 + 1 E-08x3 -1 E-05x2 + 0.0048x - 0.1933 0.1 -
R2 = 0.5226
0 -
200
400
600
800
1000
1200
1400
1600
Pressure
Fig. 8. Fraction of openable pores open vs. applied pressure for 70 x 70 grid at 40% blockage. Continuous line represents the least-squares averaged trend line.
70x70 Blockage 40% 8000
7000 6000 -
y = 0.9645x-83.534 [^ = 02412
5000 4000 3000 2000 1000 0 200
400
600
800
1000
1200
1400
1600
Pressure
Fig. 9. Flowrate vs. applied pressure for a 70 x 70 grid at 40% blockage. Continuous line represents the least-squares averaged trend line.
imum values. These values were assumed to be 1 and 10, respectively. The pseudo-random numbers are used for assigning values to the required throat diameters, so that those are uncorrelated, i.e. the size of one throat is independent of the size of any other throat. Porous media may have polymodal and/or spatially correlated pore size distributions, that can also be handled by the model, but these are not considered in this work. The flowrate vs. applied pressure relationships obtained from the model for different sets of parameters show a linear dependence similar to Darcy's law (Figs. 5,
7, 9 and 11). In these plots, the values of the y -intercept must be negative. It is indicative of the fact that only at a certain finite positive pressure does the flow through the porous medium begin to take place. The values of parameters that have been estimated using linear leastsquares technique, show a consistency in the predicted values regardless of the matrix size for a specified blockage fraction (Tables 1-5). Due to simulation time constraints we could not carry out runs at grid sizes greater than 200 for 2-dimensional networks. However looking at the graphs of the dominant para-
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
393
18x18x18 Blockage 4 0 % 1.1
pen
1.0 0.9 -
o
0.8 -
in a>
0.7 -
pen ble
o a.
0.6 -
o ro
0.4 -
fl10
0.3 -
u.
0.2 -
10
0.5 -
y= -6E-16X8 + 1E-12x5 - 1EJ)9x4 + 7E-07x3 -0.0002x2 + 0.0243x - 0.1639 R2 = 0.9738
0.1 -
nn 200
100
300
400
500
600
700
Pressure Fig. 10. Fraction of openable pores open vs. applied pressure for 18 x 18 x 18 grid at 40% blockage. Continuous line represents the least-squares averaged trend line.
18x18x18 Blockage40% 1 iUUU
10000-
8000-
£ 6000-
y=14.005x 114.95 R2 = 0.9703
4000-
2000-
0-
a
100
200
300
400
500
B00
700
Pressure Fig. 11. Flowrate vs. applied pressure for a 18 x 18 x 18 grid at 40% blockage. Continuous line represents the least-squares averaged trend line.
meter vs. grid size (Fig. 12) and cumulative flowrate vs. grid size at a specific applied pressure (Fig. 13), we believe that an infinite medium can be simulated satisfactorily in the 200-300 grid size range for 2dimensional networks. The flowrate is dependent on the number of flowpaths and that in turn on the accessible porosity of the medium. In our view, fraction of openable pores open (FOPO) is a good measure of the accessible porosity. To calculate FOPO, we noted the actual number of pores that have 'opened-up' (APO) while incrementing the pressure and then divided it by the number of pores that
are not 'blocked-off’. For 2-dimensional systems; Openable Pores = sg2(1 —fh) For 3-dimensional systems; Openable Pores = ^ ( 1 —fb) FOPO = APO=Openable Pores
(7)
(8) (9)
On plotting FOPO vs. applied pressure for different configurations, i.e. grid sizes and blockage fractions, we note an initial linear rising trend, attaining an asymptotic value subsequently (Figs. 4, 6, 8 and 10). This trend is
Table 1 Values of coefficients for fraction of openable pores open and flowrate for 2 and 3-dimensional grids at blockage of 20% Dimension
Size
Fraction of openable pores open
Flowrate R2
Coefficients
3
P5
P4
P3
P2
P1
P0
40 50 60 70 80 90 100 120 150 175 190 200
-3.00E-16 -8.00E-17 -2.00E-17 -9.00E-18 -4.00E-18 -2.00E-18 -8.00E-19 -2.00E-19 -4.00E-20 -5.00E-21 1.00E-20 1.00E-20
1.00E-12 3.OOE-13 1.00E-13 5.00E-14 2.00E-14 1.00E-14 6.00E-15 2.00E-15 5.00E-16 1.00E-16 -2.00E-16 -1.00E-16
-1.00E-09 -4.00E-10 -2.00E-10 -1.00E-10 -6.00E-11 -3.00E-11 -2.00E-11 -8.00E-12 -3.00E-12 -8.00E-13 6.00E-13 5.00E-13
7.00E-07 3.00E-07 2.00E-07 1.00E-07 7.00E-08 5.00E-08 3.00E-08 2.00E-08 7.00E-09 3.00E-09 -6.00E-10 -6.00E-10
-2.00E-04 -1.00E-04 -9.00E-05 -6.00E-05 -5.00E-05 -3.00E-05 -3.00E-05 -2.00E-05 -9.00E-06 -5.00E-06 -1.00E-06 -8.00E-07
3.28E-02 2.52E-02 2.08E-02 1.72E-02 1.46E-02 1.25E-02 1.11E-02 8.60E-02 6.40E-02 4.60E-02 2.50E-02 2.20E-02
-9.49E-02 -8.47E-01 -8.66E-01 -8.27E-01 -8.03E-01 -7.73E-01 -7.65E-01 -7.08E-01 -6.58E-01 -5.15E-01 -2.16E-01 -1.85E-01
10 11 12 13 14 15 16 17 18 19 20
-2.00E-13 -3.00E-14 -2.00E-14 -1.00E-15 -5.00E-15 -2.00E-15 -1.00E-15 -4.00E-16 -3.00E-17 1.00E-16 8.00E-17
2.00E-10 5.00E-11 3.00E-11 2.00E-11 8.00E-12 5.00E-12 3.00E-12 1.00E-12 3.OOE-13 1.00E-13 -8.00E-14
-8.00E-08 -2.00E-08 -2.00E-08 -1.00E-08 -6.00E-09 -4.00E-09 -2.00E-09 -1.00E-09 -5.00E-10 -9.00E-11 -6.00E-11
2.00E-05 6.00E-06 4.00E-06 3.00E-06 2.00E-06 1.00E-06 1.00E-06 6.00E-07 4.00E-07 2.00E-07 1.00E-07
-1.50E-03 -8.00E-04 -6.00E-04 -5.00E-04 -4.00E-04 -3.00E-04 -2.00E-04 -2.00E-04 -1.00E-04 -7.00E-05 -6.00E-05
7.16E-02 4.88E-02 4.28E-02 3.75E-02 3.14E-02 2.89E-02 2.46E-02 2.10E-02 1.77E-02 1.45E-02 1.31E-02
-2.30E-01 -1.03E-01 -8.75E-02 -8.07E-02 -2.83E-02 -4.22E-02 -2.17E-02 1.70E-03 3.88E-02 -4.50E-03 8.50E-03
Q
R2 ain el
2
P6
Coefficients P1
P0
0.9089 0.9306 0.9441 0.9539 0.9606 0.9641 0.9679 0.9746 0.9780 0.9777 0.9727 0.9748
4.89 5.98 7.81 9.28 10.51 11.65 12.88 16.79 21.10 23.58 26.98 28.20
/29.30 /40.76 /47.56 /41.79 /48.02 /81.98 /77.17 /44.28 /353.29 /259.32 /590.17 /120.83
4863.60 5935.94 7760.15 9235.21 10456.98 11567.02 12797.83 16744.72 20744.71 23324.68 26384.83 28082.17
0.8195 0.7727 0.8829 0.9014 0.9080 0.8840 0.9263 0.9195 0.9570 0.9112 0.9728 0.9568
0.9821 0.9496 0.9644 0.9680 0.9735 0.9752 0.9799 0.9805 0.9831 0.9637 0.9670
11.45 13.97 17.06 20.25 23.28 27.24 30.44 35.76 39.38 44.23 50.05
/99.17 /77.61 /96.06 /110.19 /85.03 /143.68 /164.09 /187.90 /212.43 /252.98 /283.51
11349.83 13892.40 16959.95 20136.81 23297.97 27097.32 30276.91 35573.10 39162.57 43973.02 49761.49
0.9570 0.9681 0.9729 0.9679 0.9680 0.9765 0.9786 0.9842 0.9805 0.9725 0.9684
R2 represents the regression coefficient values. Q represents the value of flowrate at an applied pressure of 1000 units. FOPO (Pa) = P6i 3 a+P5i 3 a+P4Pa + P3Pa P1Pa + P0, Q = P1 x 1000+P0.
P0, Flowrate (P a ) =
6.
!
e
a. 2
1 eg
1 ^ j
o ^~
Table 2 Values of coefficients for fraction of openable pores open and flowrate for 2- and 3-dimensional grids at blockage of 25 % Dimension
Size
Fraction of openable pores open
Flowrate
Coefficients
Coefficients
P6
P5
P4
P2
P1
P0
P1
P0
Q
R2
40 50 60 70 80 90 100 120 150 175 200
-3.00E-16 -7.00E-17 -2.00E-17 -8.00E-18 -4.00E-18 -2.00E-18 -8.00E-19 -2.00E-19 -4.00E-20 +9.00E-21 9.00E-21
8.00E-13 3.OOE-13 1.00E-13 4.00E-14 2.00E-14 1.00E-14 6.00E-15 2.00E-15 5.00E-16 -6.00E-17 -1.00E-16
-1.00E-09 -4.00E-10 -2.00E-10 -9.00E-11 -5.00E-11 -3.00E-11 -2.00E-11 -8.00E-12 -3.00E-12 -4.00E-14 4.00E-13
6.00E-07 3.00E-07 2.00E-07 1.00E-07 7.00E-08 4.00E-08 3.00E-08 2.00E-08 7.00E-09 1.00E-09 -2.00E-10
-2.00E-04 -1.00E-04 -8.00E-05 -6.00E-05 -4.00E-05 -3.00E-05 -3.00E-05 -2.00E-05 -9.00E-06 -3.00E-06 -1.00E-06
3.05E-02 2.46E-02 2.01E-02 1.70E-02 1.49E-02 1.28E-02 1.11E-02 8.60E-03 6.40E-03 3.8OE-O3 2.50E-03
-8.67E-01 -9.09E-01 -8.97E-01 -8.84E-01 + 8.95E-01 + 8.39E-01 -8.45E-01 -7.74E-01 -7.25E-01 -4.23E-01 -2.67E-01
0.8789 0.9162 0.9325 0.9444 0.9540 0.9618 0.9626 0.9710 0.9761 0.9708 0.9715
3.37 4.60 5.33 6.22 7.58 8.93 9.06 12.04 14.72 16.81 21.83
/23.16 /23.86 /34.30 /49.19 /38.44 /123.30 /54.21 /141.63 /67.14 /1.86 /103.21
3344.14 4575.55 5298.41 6174.31 7540.87 8811.10 9007.69 11902.57 14647.86 16810.14 21725.79
0.7737 0.7250 0.8331 0.7869 0.8264 0.8281 0.8748 0.9463 0.8777 0.9140 0.9471
10 11 12 13 14 15 16 17 18 19 20
-2.00E-13 -4.00E-14 -2.00E-14 -1.00E-14 -5.00E-15 -3.00E-15 -1.00E-15 -6.00E-16 -2.00E-16 1.00E-16 2.00E-16
2.00E-10 5.00E-11 3.00E-11 2.00E-11 9.00E-12 5.00E-12 3.00E-12 1.00E-12 7.00E-13 -6.00E-14 -2.00E-13
-7.00E-08 -2.00E-08 -2.00E-08 -1.00E-08 -6.00E-09 -4.00E-09 -2.00E-09 -1.00E-09 -8.00E-10 -2.00E-10 2.00E-11
1.00E-05 6.00E-06 4.00E-06 3.00E-06 2.00E-06 1.00E-06 1.00E-06 1.00E-07 5.00E-07 2.00E-07 1.00E-07
-1.50E-03 -8.00E-04 -6.00E-04 -5.00E-04 -4.00E-04 -3.00E-04 -2.00E-04 -2.00E-04 -1.00E-04 -8.00E-05 -6.00E-05
7.05E-02 5.07E-02 4.34E-02 3.83E-02 3.43E-02 2.91E-02 2.60E-02 2.23E-02 1.97E-02 1.56E-02 1.38E-02
-2.24E-01 -1.47E-01 -1.12E-01 -1.06E-01 -1.03E-01 -6.88E-02 -5.82E-02 -2.82E-02 -1.02E-02 -3.82E-02 -1.75E-02
0.9696 0.9550 0.9602 0.9673 0.9694 0.9726 0.9774 0.9813 0.9816 0.9596 0.9642
10.30 11.34 13.95 16.49 19.91 21.35 25.06 27.20 30.77 34.90 38.75
/46.54 /61.61 /79.27 /94.49 /111.98 /127.20 /146.22 /160.49 /183.64 /216.07 /244.84
10253.47 11282.40 13873.73 16392.51 19796.02 21222.80 24916.78 27038.51 30585.36 34680.93 38504.16
0.9596 0.9627 0.9470 0.9704 0.9684 0.9694 0.9800 0.9834 0.9806 0.9673 0.9660
R2 represents the regression coefficient values. Q represents the value of flowrate at an applied pressure of 1000 units. POPO (Pa) = P1Pa + P0, Q = Pl x 1000+P0.
, Flowrate (Pa) =
I3. 0
Table 3 Values of coefficients for fraction of openable pores open and flowrate for 2- and 3-dimensional grids at blockage of 30 % Size
Fraction of openable pores open
Flowrate R2
Coefficients
3
P5
P4
P3
P2
P1
P0
40 50 60 70 80 90 100 120 150 175 190 200
-2.00E-16 -5.00E-17 -9.00E-18 -5.00E-18 -2.00E-18 -9.00E-19 -4.00E-19 -1.00E-19 -2.00E-20 -1.00E-21 1.00E-20 1.00E-20
6.00E-13 2.00E-13 5.00E-14 3.00E-14 1.00E-14 7.00E-15 3.00E-15 1.00E-15 2.00E-16 6.00E-17 -1.00E-16 -1.00E-16
-7.00E-10 -3.00E-10 -9.00E-11 -6.00E-11 -4.00E-11 -2.00E-11 -1.00E-11 -6.00E-12 -2.00E-12 -6.00E-13 5.00E-13 6.00E-13
5.00E-07 2.00E-07 1.00E-07 7.00E-08 5.00E-08 3.00E-08 2.00E-08 1.00E-08 5.00E-09 2.00E-09 -5.00E-10 -9.00E-10
-2.00E-04 -1.00E-04 -5.00E-05 -4.00E-05 -3.00E-05 -3.00E-05 -2.00E-05 -1.00E-05 -8.00E-06 -5.00E-06 -1.00E-06 -2.00E-07
2.64E-02 2.14E-02 1.54E-02 1.44E-02 1.26E-02 1.09E-02 9.60E-03 8.00E-03 6.00E-03 4.40E-03 2.40E-03 1.80E-03
-7.83E-01 -8.51E-01 -7.23E-01 -8.11E-01 -7.98E-01 -8.05E-01 -7.90E-01 -7.87E-01 -7.42E-01 -6.10E-01 -2.68E-02 -1.83E-01
10 11 12 13 14 15 16 17 18 19 20
-5.00E-14 -3.00E-14 -2.00E-14 -1.00E-14 -5.00E-15 -3.00E-15 -2.00E-15 -3.00E-12 -3.00E-16 8.00E-18 1.00E-16
6.00E-11 4.00E-11 3.00E-11 2.00E-11 9.00E-12 5.00E-12 4.00E-12 6.00E-09 9.00E-13 2.00E-13 -6.00E-14
-3.00E-08 -2.00E-08 -2.00E-08 -1.00E-08 -6.00E-09 -4.00E-09 -3.00E-09 -6.00E-06 -1.00E-09 -4.00E-10 -1.00E-10
7.00E-06 6.00E-06 4.00E-06 3.00E-06 2.00E-06 2.00E-06 1.00E-06 +2.80E-03 5.00E-07 3.00E-07 2.00E-07
-8.00E-04 -8.00E-04 -6.00E-04 -5.00E-04 -4.00E-04 -3.00E-04 -3.OOE-O3 -6.85E-01 -2.00E-04 -1.00E-04 -8.00E-05
4.90E-02 5.16E-02 4.44E-02 4.01E-02 3.43E-02 3.06E-02 2.80E-02 8.28E+01 2.11E-02 1.77E-02 1.55E-02
-7.78E-02 -1.93E-01 -1.51E-01 -1.52E-01 -1.22E-01 -1.12E-01 -1.13E-01 -2.70E+02 -5.63E-01 -9.66E-02 -6.76E-02
Q
R2 ain el
2
P6
Coefficients P1
P0
0.8487 0.8784 0.8944 0.9227 0.9203 0.9383 0.9534 0.9572 0.9679 0.9728 0.9657 0.9643
2.39 3.12 3.43 4.04 4.72 5.43 6.90 7.84 10.44 13.02 12.42 13.69
/22.77 /19.96 /28.82 /34.01 /20.15 /40.26 /48.15 /39.59 /126.99 /34.94 /71.08 /144.91
2363.33 3104.14 3396.98 4004.39 4698.95 5392.14 6855.05 7804.61 10316.01 12980.06 12352.92 13544.09
0.6176 0.6818 0.7054 0.8041 0.7591 0.8179 0.8026 0.8617 0.8970 0.9100 0.8820 0.8988
0.9693 0.9671 0.9558 0.9599 0.9885 0.9732 0.9723 0.9777 0.9628 0.9581 0.9613
7.63 9.69 10.75 12.98 14.71 16.80 18.70 21.56 24.47 27.45 30.88
/68.05 /78.84 /72.81 /86.76 /98.58 /114.86 /123.03 /141.11 /137.24 /191.75 /213.58
7558.55 9615.77 10680.19 12896.24 14614.42 16684.14 18576.97 21415.89 24334.76 27262.25 30666.42
0.9357 0.9626 0.9530 0.9738 0.9621 0.9601 0.9706 0.9740 0.9765 0.9839 0.9707
R2 represents the regression coefficient values. Q represents the value of flowrate at an applied pressure of 1000 units. FOPO (PB) = P6P6a + />5i3a + P4i 3a + P3i 3 a +P2i 3 a+Pli 3 a + P0, Flowrate (Pa) = PI/* / P0, 0 = P1 x 1000+ P0
6. |
5a a.
3
1
1 •00
Dimension
^~
Table 4 Values of coefficients for Fraction of Openable Pores Open and Flowrate for 2- and 3-dimensional grids at blockage of 35 % Dimension
Size
Fraction of openable pores open
Flowrate
Coefficients
2
40 50 60 70 80 90 100 120 150 175 200
3
10 11 12 13 14 15 16 17 18 19 20
R
P6
P5
P4
P3
P2
P1
P0
-2.00E-16 -9.00E-18 -4.00E-17 -7.00E-19 -1.00E-18 2.00E-19 -4.00E-19 1.00E-19 6.00E-20 2.00E-20 1.00E-20
5.00E-13 5.00E-14 5.00E-14 7.00E-15 9.00E-15 -5.00E-16 3.00E-15 -7.00E-16 -5.00E-16 -3.00E-16 -2.00E-16
-6.00E-10 -9.00E-11 1.00E-10 -2.00E-11 -2.00E-11 -1.00E-12 -9.00E-12 1.00E-12 2.00E-12 1.00E-12 8.00E-13
3.00E-07 9.00E-08 1.00E-07 3.00E-08 3.00E-08 6.00E-09 1.00E-08 3.00E-10 -2.00E-09 -2.00E-09 -2.00E-09
-1.00E-04 -5.00E-05 -5.00E-05 -2.00E-05 -3.00E-05 -9.00E-06 -1.00E-05 -4.00E-06 -2.00E-07 1.00E-06 1.00E-06
1.95E-02 1.31E-02 1.97E-02 9.70E-03 1.00E-02 5.60E-03 6.70E-03 4.00E-03 2.20E-03 9.00E-04 9.00E-04
-5.43E-01 -4.75E-01 7.44E-02 -5.72E-01 -6.96E-01 -3.58E-01 -5.38E-01 -3.82E-01 -2.43E-01 -5.20E-03 -6.94E-02
-1.00E-13 -2.00E-14 -2.00E-14 -1.00E-14 -6.00E-15 -3.00E-15 -2.00E-15 -1.00E-15 -6.00E-16 -2.00E-16 8.00E-17
1.00E-10 3.00E-11 3.00E-11 2.00E-11 1.00E-11 6.00E-12 4.00E-12 2.00E-12 1.00E-12 6.00E-13 4.00E-15
-6.00E-08 -2.00E-08 -1.00E-08 -1.00E-08 -7.00E-09 -4.00E-09 -3.00E-09 -2.00E-09 -1.00E-09 -7.00E-10 -2.00E-10
1.00E-05 5.00E-06 4.00E-06 3.00E-06 2.00E-06 2.00E-06 1.00E-06 9.00E-07 7.00E-07 4.00E-07 2.00E-07
-1.3OE-O3 -7.00E-04 -6.00E-04 -5.00E-04 -4.00E-04 -3.00E-04 -3.00E-04 -2.00E-04 -2.00E-04 -1.00E-04 -9.00E-05
6.77E-02 4.52E-02 4.47E-02 4. HE-02 3.63E-02 3.22E-02 2.90E-02 2.57E-02 2.43E-02 2.03E-02 1.68E-02
-2.74E-01 -1.61E-01 -1.83E-01 -1.94E-01 -1.78E-01 -1.59E-01 -1.47E-01 -1.22E-01 -1.64E+01 -1.71E-01 -1.11E-01
2
Coefficients
Q
R2
P1
P0
0.6144 0.7336 0.7235 0.8182 0.8293 0.8645 0.8276 0.8983 0.9486 0.8808 0.9425
1.30 1.80 2.05 2.32 3.11 4.22 4.38 4.97 5.95 7.86 /425.30
/21.16 /16.92 /41.64 /20.71 /22.08 /19.90 /108.20 /129.83 /83.84 /483.51 8163.10
1275.94 1784.38 2010.56 2300.49 3088.02 4196.30 4268.90 4842.57 5861.36 7372.69 8163.1
0.4851 0.5058 0.5494 0.5519 0.7216 0.7665 0.6773 0.7234 0.7084 0.8720 0.8574
0.9703 0.9570 0.9499 0.9511 0.9645 0.9680 0.9690 0.9729 0.9755 0.9558 0.9565
5.62 6.63 8.18 9.54 10.97 12.78 15.49 16.74 14.01 20.71 22.93
/49.53 /41.61 /60.93 /72.61 /82.84 /92.88 /116.15 /134.24 /114.95 /162.10 /182.16
5572.47 6591.19 8120.87 9470.10 10887.16 12691.12 15377.85 16608.76 13890.05 20543.90 22747.84
0.9240 0.9507 0.9574 0.9595 0.9451 0.9546 0.9662 0.9692 0.9836 0.9772 0.9680
a.
R2 represents the regression coefficient values. Q represents the value of flowrate at an applied pressure of 1000 units. POPO (Pa) = P6Pa + P5JPa+P4JPa+P3JPa + IP2Pa 2 VlP, ,+P0, Flowrate (Pa) = PlPo / P0, 0 = P1 x 1000+ P0.
s. o o
12 ! i
3 2. eg
1
Table 5 Values of coefficients for Fraction of Openable Pores Open and Flowrate for 2- and 3-dimensional grids at blockage of 40 % Dimension
Size
Fraction of openable pores open
Flowrate
2
3
R
Coefficients P1
R2
Q
P6
P5
P4
P3
P2
P1
P0
P0
40 50 60 70 80 90 100 120 150 175
-4.00E-17 -4.00E-18 -2.00E-17 -6.00E-19 -4.00E-18 8.00E-19 -7.00E-19 3.00E-19 2.00E-19 1.00E-20
1.00E-13 6.00E-15 7.00E-14 4.00E-15 -2.00E-14 -5.00E-15 5.00E-15 -2.00E-15 -2.00E-15 -1.00E-16
-1.00E-10 9.00E-13 -1.00E-10 -1.00E-11 4.00E-11 1.00E-11 -1.00E-11 5.00E-12 6.00E-12 7.00E-13
1.00E-07 -5.00E-09 1.00E-07 1.00E-08 -4.00E-08 -1.00E-08 2.00E-08 -5.00E-09 -1.00E-08 -2.00E-09
-4.00E-05 5.00E-07 6.00E-05 1.00E-05 2.00E-05 5.00E-06 -1.00E-05 3.00E-07 1.00E-05 3.00E-06
8.40E-03 1.50E-03 1.48E-02 4.80E-03 -2.00E-03 1.00E-04 4.80E-03 2.40E-03 -6.60E-03 -7.20E-03
-7.93E-02 3.15E-01 -8.78E-01 -1.93E-01 3.63E-01 1.93E-01 -3.47E-01 -3.62E-01 1.59E+00 5.66E-01
0.4145 0.3538 0.4960 0.5226 0.5690 0.6605 0.5954 0.6712 0.6274 0.6346
0.53 1.11 0.93 0.96 1.35 2.42 1.78 3.07 3.85 5.01
/4.57 /48.56 /39.81 /83.53 /90.05 /97.95 /65.86 /871.17 /153.03 /983.12
525.53 1060.84 894.19 880.97 1262.75 2317.45 1709.44 2197.03 3697.67 4029.18
0.1420 0.4811 0.3329 0.2412 0.3091 0.3922 0.3082 0.3891 0.5664 0.5165
10 11 12 13 14 15 16 17 18 19 20
-5.00E-14 -3.00E-14 -2.00E-14 -1.00E-14 -6.00E-15 -3.00E-15 -2.00E-15 -1.00E-15 -6.00E-16 -3.00E-16 -7.00E-17
6.00E-11 5.00E-11 3.00E-11 2.00E-11 1.00E-11 6.00E-12 4.00E-12 3.00E-12 1.00E-12 8.00E-13 3.OOE-13
-3.00E-08 -2.00E-08 -2.00E-08 -1.00E-08 -7.00E-09 -5.00E-09 -3.00E-09 -2.00E-09 -1.00E-09 -9.00E-10 -5.00E-10
7.00E-06 6.00E-06 4.00E-06 3.00E-06 2.00E-06 2.00E-06 1.00E-06 1.00E-06 7.00E-07 5.00E-07 3.00E-07
-9.00E-04 -8.00E-04 -7.00E-04 -5.00E-04 -4.00E-04 -4.00E-04 -3.OOE-O3 -2.00E-04 -2.00E-04 -1.00E-04 -1.00E-04
5.27E-02 5.28E-02 4.64E-02 4.19E-02 3.78E-02 3.40E-02 3.OOE-O3 2.76E-02 2.43E-02 2.17E-02 1.91E-02
-1.81E-01 -2.68E-01 -2.34E-01 -2.45E-01 -2.35E-01 -2.27E-01 -1.96E-01 -1.94E-01 -1.64E-01 -2.19E-01 -1.85E-01
0.9378 0.9569 0.9445 0.9559 0.9605 0.9638 0.9664 0.9715 0.9738 0.9470 0.9520
3.90 5.20 5.83 7.19 8.34 9.60 11.03 12.25 14.01 16.06 17.28
/24.87 /46.21 /49.80 /59.08 /70.91 /80.99 /91.04 /101.36 /114.95 /136.46 /150.29
3870.13 5154.49 5782.70 7134.32 8266.79 9519.11 10934.97 12148.64 13890.05 15924.54 17129.71
0.9471 0.9255 0.9363 0.9373 0.9439 0.9582 0.9712 0.9659 0.9703 0.9639 0.9755
R2 represents the regression coefficient values. Q represents the value of flowrate at an applied pressure of 1000 units. POPO (Pa) = P6i3a + P5i 3 a +P4,f 1 a >P3,P a + P2 P1P a + P0, Q = Pl x 1000 + P0.
Jain
Coefficients
2
6.
, Flowrate (Pa) =
!
e
a. 2
11 . eg
1 ^ j
o ^~
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
399
Flow Rate (Slope) 100.0
4
• • •
10.0 -
•
•
•
K
•
•
n
A
A
x
•
o to
A
•
•
A
A
x
x
x
0
0
O
0
O
0
X
y.
n
x
•
• s
t
0
X
A
A X
A
A
n
A
O
0
•
o
*
1.0 -
•
•
n 20%
A
25% * 30% o 35% 140%
0.1 0
25
50
75
100
125
150
175
200
225
Grid Size Fig. 12. Slope of flowrate vs. grid size at different blockage fractions. Simulations with grid size 40 and above rate for 2-dimensional grids, while stimulations smaller than this size were for 3-dimensional grids.
FlowRate(Pa= 1000) JUUIJU
s
•
w • CgQ
10000 -
•
*
•
•
A
•A •••:
D
D
A
• A
*
•
.
n
A
y.
A
0
0
O y.
y
0
O
O
•
0
x
o
o
•
A
•
A
Q
X
1000 -
•A
• •
*
*
100 -
10 -
° 20% A 25% x 30% °35% • 40% 125
50
75
100
125
150
175
200
225
Grid Size Fig. 13. Flowrate vs. grid size at different blockage fractions at an applied pressure of 1000 units. Simulations with grid size 40 and above were for 2dimensional grids while simulations smaller than this size were for 3-dimensional grids.
expected and holds irrespective of the grid size or the fraction of pores blocked. This is similar to the observed trend when calculated cumulative mercury intrusion volume was plotted against applied pressure (Bryntesson, 2002). When the least-squares analysis is conducted and seven parameters of the polynomial estimated, it is observed that the parameter values gradually attain asymptotic values on increasing the grid size (Tables 1 -
5). This observation is in line with the fact that an ideal simulation would have all grid dimensions as infinite. At low blocked fractions, there are few chances of pores with finite resistance being surrounded by 'blocked-off pores. Thus all pores having finite resistance can be considered 'openable'. The asymptotic value of FOPO attained in all such cases is near 1 (Figs. 4 and 6). As the blockage fraction increases, the
400
S. Jain et al. / Computers and Chemical Engineering 27 (2003) 385-400
probability of finite- resistance pores being surrounded by 'blocked-off pores is higher and thus, a greater variation is seen in the asymptotic values of FOPO. This observation is most stark at blockage values within 5% of the critical value (Fig. 8). This observation can be explained by the fact that near critical values there is a sharp decrease in the probability of finding a flowpath, the implication being that there is a greater variation in the probability of finding a flowpath for a given configuration. For this to be true, it would be natural to expect a greater variation in the asymptotic value of the openable pores. These observations have been validated by the regression analysis while computing the parameter values. The R2 trend values indicate that there is greater deviation from 1 at blockage fractions near critical values, particularly at low grid sizes (Tables 1-5). As grid size increases the R2 values improve indicating that the simulation is closer to that of an infinite medium.
5. Conclusions
This simulation is a simplified but innovative way of dealing with flow through porous media. Using the Monte Carlo technique and concepts of percolation theory, simulations have been carried out for 2 and 3dimensional models of porous media in capillary flow regime. The flowrate vs. applied pressure results have been found to be in agreement with Darcy's law. Parameter estimation of fraction of openable pores that have 'opened up' vs. applied pressure has also been shown to be consistent. More accurate results could be obtained by using actual resistance distributions of porous media. It might then be possible to predict flowrate vs. applied pressure behaviour for systems with different physical properties (e.g. surface tension and contact angle). Actual experimental data for flow in the capillary regime have not been reported in literature, where the major focus has been on statistical description of porous media. Also, experiments that reproduce the exact conditions of the simulation would be difficult to perform. Therefore, the results mentioned herein should be treated as being of a qualitative nature. The parameter values have been computed after considering several thousands of points for each configuration (matrix size and blocked fraction). Runs were concluded once the curve predicting FOPO vs. applied pressure exhibited an asymptotic behaviour. It has been our observation that until this asymptotic behaviour is observed the runs should not be concluded, otherwise the parameter values would be misleading. This is also in accordance with the random nature of the simulation
because theoretically an infinite number of runs are required before presenting any trend. Further work in this field would include extension of the simulation to the region of Hagen-Poiseulle flow, where the resistance is inversely related to the fourth power of the radius of the element. The resistance beyond the laminar flow region decreases slightly and then reaches a constant value at very high flowrates in the turbulent zone. This can be modelled as two separate flow regions, with elements having different (but constant) resistances in each region.
Acknowledgements
We are grateful to Dr M.N. Gupta, Senior Manager, Computer Services Centre, IIT, Delhi, for the generous allocation of computational time. We would also like to thank all staff members of CSC, IIT, Delhi, for their cooperation throughout this project. In particular, we would like to mention Mr. Gopal Krishen, the system programmer, and Mr. Gulshan Naveriya, the hardware engineer, for their invaluable assistance.
References Bernsdorf, J., Brenner, G., & Durst, F. (2000). Numerical analysis of the pressure drop in porous media flow with lattice Boltzmann (BGK) automata. Computer Physics Communications 129, 247255. Binder, K. (Ed.), Monte Carlo methods in statistical physics (1979). Springer-Verlag. Bryntesson, L. M. (2002). Pore network modelling of the behaviour of a solute in chromatography media: transient and steady-state diffusion properties. Journal of Chromatography A 945, 103-115. Cordero, S., Rojas, F., & Riccardo, J. L. (2001). Simulation of threedimensional porous networks. Colloids and Surfaces A: Physicochemical and Engineering Aspects 187-188, 425-438. Cormen, T. H., Leiserson, C. E., & Rivest, R. L. (2000). Introduction to algorithms. EEE: Prentice Hall of India. De Weist, R. J. M. (Ed.), Flow through porous media (1969). Academic Press. Dullien, F. A. L. (1992). Porous media—fluid transport and pore structure. Academic Press. Efros, A. L. (1982). Physics and geometry of disorder—percolation theory. Mir Publishers. Greenkorn, R. A. (1983). Flow phenomenon in porous media. Marcel Dekker. Hilpert, M., & Miller, C. T. (2001). Pore-morphology-based simulation of drainage in totally wetting porous media. Advances in Water Resources 24, 243-255. Manegold, V. E. (1937). Uber Kapillarsysteme. Kolloid-Z 80, 253. Monteagudo, J. E. P., Rajagopal, K., & Lage, P. L. C. (2002). Simulating oil flow in porous media under asphaltene deposition. Chemical Engineering Science 57, 323-337. Muskat, M. (1982). Flow of homogenous fluids through porous media.
IHRDC. Scheidegger, A. E. (1963). The physics of flow through porous media. University of Toronto Press.
