Two-phase fluid flow in geometric packing.

Downloaded from http://rsta.royalsocietypublishing.org/ on April 10, 2016

Two-phase fluid flow in geometric packing rsta.royalsocietypublishing.org

Aureliano Sancho S. Paiva1 , Rafael S. Oliveira1,2 and Roberto F. S. Andrade1

Research Cite this article: Paiva ASS, Oliveira RS, Andrade RFS. 2015 Two-phase fluid flow in geometric packing. Phil. Trans. R. Soc. A 373: 20150111. http://dx.doi.org/10.1098/rsta.2015.0111 Accepted: 27 July 2015 One contribution of 13 to a theme issue ‘Topics on non-equilibrium statistical mechanics and nonlinear physics (II)’. Subject Areas: fluid mechanics, complexity Keywords: bi-phase flow, Apollonian packing, oil drainage Author for correspondence: Roberto F. S. Andrade e-mail: [email protected]

1 Instituto de Física, Universidade Federal da Bahia, 40210-210

Salvador, Brazil 2 Centro de Formação de Professores, Universidade Federal do Recôncavo da Bahia, 45300-000 Amargosa, Brazil We investigate how a plug of obstacles inside a two-dimensional channel affects the drainage of high viscous fluid (oil) when the channel is invaded by a less viscous fluid (water). The plug consists of an Apollonian packing with, at most, 17 circles of different sizes, which is intended to model an inhomogeneous porous region. The work aims to quantify the amount of retained oil in the region where the flow is influenced by the packing. The investigation, carried out with the help of the computational fluid dynamics package ANSYSFLUENT, is based on the integration of the complete set of equations of motion. The study considers the effect of both the injection speed and the number and size of obstacles, which directly affects the porosity of the system. The results indicate a complex dependence in the fraction of retained oil on the velocity and geometric parameters. The regions where the oil remains trapped is very sensitive to the number of circles and their size, which influence in different ways the porosity of the system. Nevertheless, at low values of Reynolds and capillary numbers Re < 4 and nc 10−5 , the overall expected result that the volume fraction of oil retained decreases with increasing porosity is recovered. A direct relationship between the injection speed and the fraction of oil is also obtained.

1. Introduction Flows of bi-phase fluids in porous media occur in a large number of industrial applications and in natural systems. Among them, we may indicate boiling and pressurized water reactors, phase separators used in oil production, geo-meteorological phenomena, including sedimentation and soil erosion, blood flow [1–4]. 2015 The Author(s) Published by the Royal Society. All rights reserved.


We consider a two-dimensional channel in the plane XY, limited by two parallel lines placed at y = ±h/2, where the fluid flow takes places under a pressure gradient applied along the

.........................................................

2. Model

2

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 373: 20150111

As well known, the phenomenological concept ‘porous medium’ encompasses a multitude of physical conditions that share a common property: the presence of solid matter within the region where the flow is observed. The medium can be in the form of a single plug, inside which many narrow channels allow for the fluid flow, or in the form of a large number of solid disconnected obstacles, which also lead to the presence of very narrow channels as compared to the linear dimensions of region [5]. The first situation often occurs in the subsurface, where large pressures compact geological layers to form porous sedimentary rocks, which form the preferential storage system for petroleum, gas and water. The second situation can be thought of as sand and gravel systems on river bed, but also on subsurface systems where the sedimentary rocks have been fractured either by natural or anthropogenic factors [6,7]. Both situations were first described in a uniform way, by means of the permeability k, a single phenomenological parameter characterizing the global effect of the solid matter that forces the fluid to flow along a much more complex pattern, with much higher energy dissipation when compared with the corresponding obstacle-free flow. Here k depends on several geometrical and physical properties of the medium, like the porosity φ [3,4]. The investigation of the dependence of k on the actual distribution of solid matter largely benefited from computational fluid dynamics (CFD) methods [8,9], i.e. from the direct integration of flow equations. Another strategy in the study of fluid motion that also received much attention consists in replacing the actual medium by cellular models [10]. Macroscopic quantities are treated locally, by assuming that the displacement of the fluid among neighbouring cells must obey the proper conservation laws. In the case of porous media, this amounts to consider a discrete assembly of interconnected channels, and assuming that within each channel the flow is obstacle free [11,12]. Both possibilities, initially proposed to single-phase problems, can be conveniently extended to consider multi-phase flows. Here, the scope of the problem is enhanced, and questions beyond understanding the details of resistance of the medium to the flow arise. The most direct one is to describe how the obstacles change the regions that different fluids occupy. In this work, we consider the flow of a two-phase Newtonian fluid (e.g. water and oil) in a two-dimensional channel that accommodates a number of circular obstacles (discs), geometrically displaced according to the Apollonian packing [13–16]. This particular assembly, which has been used in previous investigations [17,18] about the dependence of k with respect to packing composed of obstacles of different sizes, aims to describe the main properties of a fractured medium with grains of different sizes. For such packing, φ and k depend on the number and size of discs [17,18], a fact that influences the results discussed herein. We would like to call attention that usual porous medium studies consider that the whole region where the fluid flows is somewhat homogeneous. This differs from the described model set-up, which is characterized by the presence of obstacles in a limited region of the channel, the rest being obstacle free. This limitation is due to our intention to concentrate the study to understanding the effect of a number of obstacles in retaining some fraction of a resident (defending) viscous fluid (oil) when the region is invaded by a less viscous fluid. The main goal of this study is to understand how the presence of the packing limits the access of each phase to different regions of inside the channel. Since we consider, as initial conditions, that oil occupies the space where the packing is placed, we obtain information about the fraction of oil retained inside the channel due to the presence of the packing. Our investigation explored the dependence of the retained oil with respect to the input velocity and the two geometric aspects (size and number of discs) defining the porosity of the packing. Since no analytic solution of the Navier–Stokes for this problem is available, it can only be addressed within CFD or other numerical approach as, for instance, cellular automata and lattice Boltzmann methods.


(a)

3

L 3

4

h

(b)

Figure 1. (a) First and (b) third generations of the Apollonian packing inside the channel of height h and extension L. The numerical integration considers non-slip boundary conditions in the interface between any solid object (all discs in the packing and channel walls) and any fluid. (Online version in colour.)

X-direction. In its interior, a set of discs, grouped according to an Apollonian packing, models a fractured porous region, composed of isolated obstacles. This packing was introduced by Apollonius of Perga in the third century BC as the solution to the problem of determining the minimum number of non-overlapping circles necessary to cover a region in the plane [19]. In our model, the packing is constructed according to the following rule. Initially, at the first packing generation G = 1, four tangent discs of the same size, with radius R = h/4, are placed together within the channel, in such a way that their centres are located on the vertices of a square of side h/2. The second generation G = 2 is obtained by adding three more discs: one in the centre of the G = 1 packing, which touches the four G = 1 discs, and two discs in the space limited by two discs and one of the walls, which are considered to be circles with infinite radii placed at y = ±∞. Higher order generations G > 2 are obtained by adding up new tangent circles in all spaces limited by any three discs in the previous generation. The number of discs in the packing at generation G is indicated by NG . We recall that, given any set of three mutually tangent circles with known radii and centre locations, it is possible to determine the radii and the coordinates of the centres of a fourth circle that is tangent to the three previous ones [19]. The touching discs in the described configuration refrain fluid movement under the pressure gradient for any value of G. To allow for fluid displacement, it becomes necessary to open channels between the discs, which is achieved by uniformly reducing their radii Ri by a factor s < 1. In this work, we consider s = 0.7, 0.8 and 0.9, and G = 1, 2 and 3. In this way, the packing porosity, which takes account of the space between the discs within a square of area h2 , depends both on s and on G, i.e. φ = φ(s, G). For the purpose of finding the numerical solution within the CFD procedure, the otherwise infinite channel is considered to have a finite length L. The model set-up is illustrated in figure 1, for G = 1 (a) and G = 3 (b), and s = 0.9. For the purpose of keeping track of the presence of the defending (oil) and invading (water) fluids in different regions of the channel, it is divided into four zones Zi i = 1 4, indicated by different colours (grey tones). With this division, it is possible to use a facility in the CFD application to choose the area in which we want to calculate the volume of each of the phases. At the initial integration time (t = t0 = 0), the invading fluid (water) occupies zone Z1 . The void space of zone Z2 , as well as zones Z3 and Z4 are fully filled with oil. This model has already been used to investigate the flow of single-phase fluid under different boundary conditions. Although the particular choice of disc sizes and geometrical placement can hardly be found in natural systems, the use of the Apollonian packing to model porous geometry has some advantages, which result from this specific packing geometry. Indeed, among many other properties, it is observed that the disc sizes decay exponentially as their number in the packing grow, the distribution of discs with respect to their number of neighbours in the packing decays as a power law, the network obtained connecting the centres of touching discs, which reproduces the structure of contact forces, is topologically invariant. When it is used in any numerical modelling, the packing is limited to a finite number of sites but, due to

.........................................................

2


1


The equations of mass and momentum conservation are integrated with the help of ANSYSFLUENT CFD package, under the assumption of incompressible and isothermal regime [22]. The equation of conservation of mass, or the continuity equation, is written as ∂ρ + ∇ · (ρv) = Sm . ∂t

(3.1)

This equation is valid for both compressible fluids as for incompressible fluids, and the term Sm is the flow mass added to the continuous phase from a dispersed phase. The incompressibility assumption requires that ρ does not depend on time. Thus, if Sm = 0, the continuity equation reduces to ∇ · v = 0. (3.2) The equation for momentum conservation, in an inertial system of reference, is written as ∂ (ρv) + ∇ · (ρvv) = −∇p + ∇.(τ¯ ) + ρg + F, ∂t

(3.3)

where p is the pressure, τ¯ is the stress tensor, ρg is the gravitational force on the fluid mass and F represents any external force acts in the fluid. The stress tensor τ¯ is given by τ¯ = μ[(∇v + ∇vT ) − 23 ∇ · vI],

(3.4)

where μ is the molecular viscosity, I is the unit tensor and the second term on the r.h.s. describes the effect of volume expansion. Thus, ignoring external forces, and using the non-slip condition and the continuity equation in equations (3.3) and (3.4), the momentum conservation equation is reduced to the Navier–Stokes equation ρ

∂v + ρv. ∇v = −∇P + μ∇ 2 v. ∂t

(3.5)

The numerical integration takes into account the usual non-slip boundary condition (v = 0) between all solid–fluid interfaces. The CFD framework requires the definition of a mesh, which

.........................................................

3. Methods

4


its deterministic nature, it is easy to improve the quality of the results, just by increasing the value of G (and number of discs) in a controlled way. As we mentioned in the Introduction, the geometrical properties of the packing provide an interesting way to analyse the properties of flows in a region occupied by objects of different sizes, for instance, a fractured rock. The model can also be changed by randomly disturbing the geometrical setting. This amounts to changing the positions and/or sizes of discs, as well as individually selecting the reduction factor for each object. By controlling amplitude of such changes, it is possible gradually approach to more realistic description of porous media formed by isolated objects [20]. Finally, the number of discs can also be increased by adding more than one set of discs defined by the Apollonian packing, which can be subjected to random disturbances as discussed earlier. Of course, the investigation of all possible generalizations is far beyond the scope of this work. Therefore, we concentrate on the results for the pure deterministic packing. In such case, it is not necessary to take averages over discs distributions. Nevertheless, we have performed several independent integrations for the same geometry, as the CFD results may be susceptible to uncontrolled rounding-off error propagations. The geometrical set-up considered the channel width h = 4.0 mm, while the total length L = 33 mm is assigned to the four zones according to the following scheme: Z1 : 1 mm; Z2 : 4 mm (packing zone); Z3 : 8 mm; Z4 : 20 mm. By properly re-scaling the linear dimensions, velocities and fluid features of the fluid used in this study, our results can also be applied to similar systems of different sizes. The fact that the permeability of similar systems of different size does not differ much has been pointed out as an advantage of numerical methods in other approaches to porous media problems [20,21].


vx t vy t + ≤ Cmax , x y

(3.6)

where C is the (dimensionless) Courant number, t is the time step, vx (vy ) is the velocity along the x (y)-direction, and x (y) is the displacement of the fluid phase travelled within the cell along x (y)-direction. The CFL condition relates the typical time interval required for the fluid parcel to traverse any given cell. If the CFL condition is not satisfied, errors associated with unphysical deformation of the two fluid interface may propagate, so that the numerical solution becomes meaningless. For two-phase systems, condition (3.6) must be satisfied for the velocities of both fluids, which leads to difficulties in adjusting space discretization and time step in order to optimize the performance of the package. Indeed, when ANSYS-FLUENT performs a time-dependent multi-phase fluid flow calculation [22], the time step used for the volume fraction calculation will not be the same as the time step used for the rest of the transport equations. In the region near the fluid interface, the software divides the area of each cell by the sum of the outgoing fluxes. The result represents the time it would take for the fluid to empty out that cell (with known size). The characteristic transit time is taken as the minimum value of this time taken over all the cells. It is used to verify whether condition (3.6) holds.

4. Results The results we obtained were limited to the low Reynolds number regime (Re < 10), where Re =

hvρ ¯ . μ

(4.1)

For this range of Re values, the properties of the flow also depend on the capillary number nc , defined by μv¯ . (4.2) nc = σ In equations (4.1) and (4.2), h indicates the channel height, v¯ the constant inlet speed of the invading fluid, ρ and μ represent the invading fluid density and viscosity, respectively, while σ is the surface tension between the two fluids [5,23]. In this investigation, the capillary number, which expresses the relationship between the viscous force and surface tension between the two fluids, takes values in the interval 3.0 × 10−6 < nc < 3.5 × 10−5 . So it covers not only a range conditions where capillary forces dominate viscous forces nc ∼ = 10−5 , but also the opposite situation in which viscous forces dominate the flow.

.........................................................

C=

5


results from a division of the domain into discrete, non-overlapping cells [9]. The differential equations expressing mass and momentum conservation are transformed into algebraic equations for the values of the pertinent fields on the faces of the cells. The resulting values for each of the two fluid phases are interpolated to the centre of cells. For a model formed by circles, as the Apollonian packing, it is advisable to use a discretization based on triangular cells, which allow for a best fit of the curved walls. The mesh generator used to discretize the domain ignores the symmetry of the packing geometry, resulting in a triangulation, which does not obey mirror symmetry with respect to the centre line. In this study, the end effect is that the mesh does not exactly reproduce the reflection symmetry with respect to the horizontal line at centre of the channel. In the next section, it is possible that asymmetric oil retention patterns result as stationary state at the end of invading process. Therefore, a completely symmetric mesh would be relevant in this study in order to discern if any pattern asymmetry results from discretization or from the intrinsic instability of a more symmetric solution. Nevertheless, we were able to greatly circumvent the mesh asymmetry by imposing the same mesh boundaries at all fixed walls, which significantly reduced mesh distortions. In the numerical solution of equations (3.2) and (3.5) for two-phase systems, the Courant– Friedrichs–Lewy (CFL) condition is taken in account to analyse the stability of the solution [9] and controlling its accuracy. For two-dimensional systems, the CFL condition is expressed by


Table 1. Velocity, Reynolds number and capillary number. −1

6

nc

0.0001

0.39

3.21 × 10−6

0.00025

0.99

8.03 × 10−6

0.0005

1.99

1.60 × 10−5

0.00075

2.98

2.41 × 10−5

0.001

3.98

3.21 × 10−5

.......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... ..........................................................................................................................................................................................................

The effect of the packing on the drainage of oil from the channel by the less viscous fluid depends on the values of the following parameters: contact angle, inlet speed, channel size, fluid viscosities and surface tension [23]. For each of the quoted values of reduction factor s and generation G, we carried on numerical integration of the system for five different values of inlet speed v¯ indicated in table 1. During this procedure, the values of all other parameters were kept fixed. For the invading fluid, we used the properties of water at 20◦ C, which corresponds to ρinv = 998.2 kg m−3 and μinv = 0.001003 kg ms−1 , while the defending fluid is characterized by ρdef = 922 kg m−3 and μdef = 0.05 kg ms−1 . This corresponds to oil with typical density (in respect to that of water) expressed by American Petroleum Institute (API) index 22. Further parameter values are σ = 0.0312 N m−1 , contact angle between fluids θ = π/2, and time integration step τ = 2.5 × 10−4 s. The choice of the inlet speed range was based on the analysis of the flow of the same two-phase system in a channel without any obstacle. We selected values of v¯ such that, in the packing-free condition, all oil could be expelled by the invading water. In this way, we can be sure that any oil retention was caused only by the presence of the packing. The corresponding value of Re were low enough to be in the laminar (table 1). There is an intermediate range of values of the inlet velocity in which the regime is still laminar but the invading fluid is not able to completely expel the oil from the channel. This regime is characterized by a particular fingering pattern, in which a single finger is opened connecting inlet to outlet, through which water flows. In the stationary regime, a certain amount of oil remains adhered to the channel walls [23]. We measure the effectiveness of the oil drainage by the quantities B2 and B2,3 . Let Oi , i = 2, 3, denote the amount of oil retained in zones Z2 and Z3 . The amount of oil initially placed in the G 2 zone Z2 corresponds to A2 = h2 − N i=1 π ri (s), where ri (s) indicates the radius of each disc in the packing with reduction factor s. Then let us define B2 =

O2 A2

and B2,3 =

O2 + O3 , A2

(4.3)

so that B2 informs how much oil was expelled from Z2 by the invading fluid, and B2,3 measures the total amount of oil that could not be recovered at the end of the process.

(a) Process of oil retention When water is injected at a constant speed, it first pressures the amount of water initially placed in Z1 against oil in Z2 , forcing it into Z3 . After some time interval, water also reaches Z3 but leaves behind a certain amount of oil between the discs in the zone Z2 . Depending on the features of the two fluids and on the value of v, ¯ the channels opened by the invading fluid (consequently, the regions where oil remains stored) may have different patterns. The fraction of retained oil depends essentially on the resulting channel pattern. Variations on the retained fraction (e.g. when v¯ is increased) are much smaller when the channel pattern is not changed when compared with the opposite condition, in which a new channel pattern arises. In all investigated situations, water

.........................................................

Re


Vel (m s )


(a)

Figure 2. Typical trapped oil pattern in the first-generation packing for (a) v¯ = 0.001 m s−1 , (b) v¯ = 0.00075 m s−1 and (c) v¯ = 0.0005 m s−1 . Patterns in (b,c) do not display inversion symmetry. (Online version in colour.)

flow through the central channel, horizontally placed between the large G = 1 discs, has been observed. The central channel is the preferential path taken by the invading fluid. Water reaches Z3 in three different patterns: (i) exclusively through the central channel; (ii) through the central channel and one of the two channels formed by a large circle and the wall; and (iii) through all three channels. In this last case, the flow through Z3 and Z4 is much like that observed in the packing-free channel at the same value of v. ¯ The invading fluid completely expels oil from these regions. The amount of retained oil (if any) is trapped within the channels between the discs and, eventually, a small drop attached to the right part of the discs in the downstream direction. In case (i), the invading fluid flowing through the central channel pressures oil towards the walls. For the chosen values of v, ¯ the invading fluid causes the progressive decrease of the width of the oil layers adhered to the walls, until the rupture of the resulting thin films occurs. This happens at a certain distances (not necessarily the same for the two walls) from the point marking the beginning of Z3 . At the moment this occurs, the large surface tension causes the oil left behind to flow back in the −x-direction. An amount of oil remains attached to the large circle, part in zone Z2 , part in Z3 . Similarly, the region occupied by the downstream oil is also reshaped towards the usual interface pattern of the packing-free invasion process. Finally, the dynamics in case (ii) follows a combination of the processes observed in (i) and (iii). The different scenarios described above depend on the number and size of discs through changes in the total area of the walls of the channel, which impact on the influence of the viscous forces. It is possible to identify a trend towards the enhancement in the amount of oil expelled with decreasing speed, which results from changes between capillary forces and viscous forces. The changes in oil retention profile are also related to accumulation in the downstream region after the packing, as can be seen in figures 2 and 3. The qualitative analysis of the flow by the direct comparison of resulting patterns indicates that the presence of new circles causes the water to penetrate the porous region, as can be seen by the comparison between figures 2 and 3. These illustrations indicate that, in a laminar flow with low Re, despite the larger energy dissipation, owing to the increase in the fluid–wall interface caused the increase of the number of discs, geometric aspects of the porous region can facilitate the invasion of water into thinner channels. However, this is not a guarantee of an enhancement of oil removal, because oil can be trapped in other regions. These can be better understood by combining the observation of trapping patterns with quantitative measures expressed by the parameters B2 and B2,3 .

.........................................................

(c)


(b)

7


(a)

8

Figure 3. Typical trapped oil patterns for G = 2, s = 0.8 and (a) v¯ = 0.0005 m s−1 , (b) v¯ = 0.00075 m s−1 and (c) v¯ = 0.001 m s−1 . All patterns lack inversion symmetry with respect to the centre of the channel. (Online version in colour.) (b)

60

60

50

50

40

40

30

30

20 60

20 60

50

50 B2,3

B2

(a)

40

40

30

30

20 60

20 60

50

50

40

40

30

30

20

20 0

0.02

0.04 0.06 0.08 v–(cm s–1)

0.10

0

0.02

0.04 0.06 0.08 v–(cm s–1)

0.10

Figure 4. Volume fractions (a) B2 and (b) B2,3 as function of v. ¯ Squares, circles and triangles correspond, respectively, to G = 1, 2 and 3. From top to bottom, s increases from 0.7 to 0.8 and 0.9. (Online version in colour.)

(b) Analysis of the volume fraction of retained oil Figure 4 illustrates the dependence of B2 and B2,3 with respect to v¯ for G = 1, 2, and 3 and s = 0.7, 0.8 and 0.9. It shows that B2 remains almost constant for different values of v¯ when G = 1. This is consistent with the illustrations in figure 3, showing that the invading fluid flows predominantly through the central channel, while the four secondary channels between circles

.........................................................

(c)


(b)


60

(b)

60 40

30

30

20 60

20 60

50

50

40

40

30

30

20 60

20 60

50

50

40

B2,3

50

40

40

30

30

20 60

20 60

50

50

40

40

30

30

20 60

20 60

50

50

40

40

30

30

20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 f(s, g)

20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 f(s, g)

Figure 5. Volume fractions (a) B2 and (b) B2,3 as function of φ(s, G). Squares, circles and triangles correspond, respectively, to G = 1, 2 and 3. For fixed G, φ decreases for increasing values of s = 0.7, 0.8 and 0.9. From top to bottom, v¯ assumes the following values (in cm s−1 ): 0.01, 0.025, 0.05, 0.075 and 0.1. (Online version in colour.)

and walls remain blocked. The only exception observed for the smallest value of v¯ when s = 0.7 results from opening one channel between one of the left discs and the wall. On the other hand, for the same values of s and v, ¯ the dependence of B2,3 is subject to larger variations. This results from the different amounts oil that remain attached to the right of the discs, most of which occupy regions in zone Z3 . This picture changes for the higher generations, mainly because of the presence of the disc with the second largest radius in the central part of the packing. This disc alters the flow pattern in different ways for different values of s and v. ¯ The dependence of B2 as function of v¯ is characterized by large variations for all values of s, some of which with opposite tendencies. Fluctuations in the behaviour of B2,3 are enhanced, since it results from oil stacked within the small channels as well as in the region behind the large circles zone Z3 . ¯ an overall general Despite the large fluctuations in the values of B2 and B2,3 with respect to v, tendency is the increase in the fraction of retained oil as the value of s increases. Although this aspect becomes more evident for B2,3 , it can also be detected in the case of B2 . This corresponds to the reported behaviour in other previous studies and the overall known resistance to remove a fluid from dense packing and media with small porosity. The larger differences between B2 and B2,3 when the value of s increases is explained by the fact that the amount of oil O3 is less sensitive to variations in the value of s than O2 and A2 . It is important to note that the two different measures can guide expected results for other packing arrangements that can be built in similar way. If the whole channel is filled with discs, one should expect the value of B2 as a good approximation. However, for inhomogeneous packing, where regions with obstacles and obstacle-free alternate, B2,3 should provide a better

.........................................................

50

9


B2

(a)


The results for the fraction of retained oil can also be analysed as a function of the packing porosity, i.e. as function of φ = A2 /h2 . This is illustrated in figure 5, where we show graphs for ¯ Once φ = φ(s, G), each panel B2 and B2,3 as a function of φ for each of the five different values of v. comprise nine points. For the sake of a better visualization, we identify the points corresponding to the same value of G by different symbols and/or colours, connected by straight lines, so that the variation of φ results from different values of s. The graphs for B2,3 show a decrease in the fraction of retained oil with increasing porosity, that is, it becomes easier to expel the oil from the packing for larger values of φ, when the channels between the discs and between the discs and the walls are wider. This dependency becomes weaker in the case we compare the behaviour of B2 . For instance, when G = 1, the fraction of oil retained within the packing region can be regarded as rather insensitive to the value of s for almost all values of v. ¯ For G = 2 and 3, the tendency of a decrease in B2 with respect to φ is again recognized, but the range of retained values is clearly smaller. As already discussed, the changes in the trend of the reduction of retained oil with respect to s are related to different accumulation patterns indicated in figures 2 and 3. This occurs more often at low values of v. ¯ The graphs also make evident that, at G = 1 there is less variation, what can be explained by the dominating role of the pattern where only the central channel is opened by the invading fluid.

5. Conclusion This work discussed the numerical investigation of a two-phase flow. The main focus was to obtain a quantitative relationship between the amount of oil retained within a two-dimensional channel with a packing of discs of different sizes in its interior when the system is invaded by a less viscous fluid. Although the flow differs from that occurring in a Hele–Shaw set-up by neglecting the viscosity effect between fluids and limiting plates, viscous fingering phenomena are also present [24]. The packing is constructed by a geometric recurrence rule, and the flow proceeds at low speed to warrant a laminar flow. The regime of low velocity was also reduced to warrant that, at all invading velocities used in the investigation, no oil would be retained inside the channel in the packing-free condition. The total amount of retained oil depends on the invading fluid velocity as well as on several geometrical features, including the number, position and size of the discs. No simple relation could be obtained between this amount and any single specific geometric feature of the packing. Global tendencies, however, can be detected. They are consistent with empirical findings that extraction is difficult when porosity is small and the invading velocity increases. The quantitative results were expressed in terms of the fraction of oil retained in the region where the packing was contained, as well as in terms of all remaining oil. This makes it possible to extrapolate the results to situations in which a larger number of obstacles is uniformly placed within the channel, or when they follow a sequence of regions with and without obstacles. Different effects follow the increase of the number of discs inside the channel. On one side, it creates barriers within the porous region and increases the area of contact of fluid with the walls. At the same time, changes in the flow pattern may enhance the efficiency of the invading fluid in expelling the resident one. Preliminary results for more complex configurations, with two or more Apollonian plugs placed inside the channel, confirm that the devised framework can lead to similar results for more realistic patterns (see §2), which will be discussed in a forthcoming contribution. The computational efforts to investigate many distinct configurations are much larger, as we have

.........................................................

(c) Porosity

10


approximation for the fraction of retained oil. This is useful when one compares the dependency of these measures as function of v¯ as well as for the analysis of the dependency with respect to the porosity, that is carried out in the following subsection.


References 1. Scheidegger AE. 1974 The physics of flow through porous media, 3rd edn. Toronto, Canada: University Toronto Press. 2. Bear J. 2013 Dynamics of fluids in porous media. New York, NY: Dover Civil and Mechanical Engineering, Dover Publications. 3. Sahimi M. 2011 Flow and transport in porous media and fractured rock: from classical methods to modern approaches, 2nd edn. Weinheim, Germany: Wiley-VCH. 4. Dullien FAL, Brenner H. 2012 Porous media: fluid transport and pore structure. Amsterdam, The Netherlands: Elsevier Science. 5. Tiab D, Donaldson EC. 2012 Petrophysics: theory and practice of measuring reservoir rock and fluid transport properties, 2nd edn. Burlington, Gulf Professional Publication. 6. Vandenbossche J, Yox-theimer D, Vidic R, Brantley S, Abad J. 2013 Impact of shale gas development on regional water quality. Science 340, 1235009. (doi:10.1126/ science.1235009) 7. Male F, Patzek T, Marder M. 2013 Gas production in the Barnett Shale obeys a simple scaling theory. Proc. Natl Acad. Sci. USA 110, 19 731–19 736. (doi:10.1073/pnas.1313380110) 8. Versteeg HK, Malalasekera W. 2007 An introduction to computational fluid dynamics: the finite volume method. Harlow, UK: Longman Scientific and Technical. 9. Tu J, Yeoh GH, Liu C. 2012 Computational fluid dynamics: a practical approach. London, UK: Butterworth-Heinemann. 10. Rothman DH, Zaleski. S. 1994 Lattice-gas models of phase separation: interfaces, phase transitions, and multiphase flow. Rev. Mod. Phys. 66, 1417–1479. (doi:10.1103/RevModPhys. 66.1417) 11. Wilkinson D, Willemsen JF. 1983 Invasion percolation: a new form of percolation theory. J. Phys. A 16, 3365–3376. (doi:10.1088/0305-4470/16/14/028) 12. Ferreol B, Rothman DH. 1995 Lattice-Boltzmann simulations of flow through fontainebleau sandstone. In Multiphase flow in porous media (ed. PM Adler), pp. 3–20. Berlin, Germany: Springer. 13. Adler P. 1985 Transport processes in fractals. I-Conductivity and permeability of a Leibniz packing in the lubrication limit. Int. J. Multiphase Flow 11, 91–108. (doi:10.1016/03019322(85)90007-2) 14. Adler PM, Thorvert JF. 1998 Real porous media: local geometry and macroscopic properties. Appl. Mech. Rev. 51, 537–585. (doi:10.1115/1.3099022) 15. Doye JPK, Massen CP. 2005 Self-similar disk packings as model spatial scale-free networks. Phys. Rev. E 71, 016128. (doi:10.1103/PhysRevE.71.016128) 16. JS Andrade, Herrmann HJ, Andrade RFS, da Silva LR. 2005 Apollonian networks: simultaneously scale-free, small world, euclidean, space filling, and with matching graphs. Phys. Rev. Lett. 94, 018702. (doi:10.1103/PhysRevLett.94.018702)

.........................................................

Authors’ contributions. A.S.S.P. conceived the study, prepared the CFD framework, carried out numerical work and result analyses, and helped draft the manuscript; R.S.O. conceived the study and prepared the CFD framework; R.F.S.A. conceived the study, coordinated the study and wrote the manuscript. All authors gave final approval for publication. Competing interests. We declare we have no competing interests. Funding. The authors acknowledge the financial support of Brazilian agencies CNPq, CAPES, FINEP, and FAPESB. R.F.S. Andrade benefit from the support of the Instituto Nacional de Ciência e Tecnologia para Sistemas Complexos (INCT-SC).

11


to consider all difficulties that this study made evident for the simulation of two-phase flows within CFD. A good mesh discretization becomes crucial due to the necessity to account for the fractions of two fluids in the same cell. The control of Courant number, which is not relevant to the single-phase flow in the same packing, becomes an important issue. For a fixed mesh, this control refrains the range of velocities of the invading fluid for which the condition expressed by equation (3.6) is not violated. Finer meshes would be required to reach larger velocities, which constitute a limitation to the investigation as compared to the single-phase system.


12 .........................................................


17. Oliveira RS, Andrade JS, Andrade. RFS. 2010 Fluid flow through Apollonian packings. Phys. Rev. E 81, 047302. (doi:10.1103/PhysRevE.81.047302) 18. Oliveira RS, Andrade JS, Andrade. RFS. 2015 Fluid flow through packings of rotating obstacles. Phys. Rev. E 91, 033002. (doi:10.1103/PhysRevE.91.033002) 19. Mallows CL, Lagarias JC, Wilks AR. 2002 Beyond the descartes. Circle theorem. Am. Math. Mon. 109, 338–361. (doi:10.2307/2695498) 20. Araújo AD, Bastos WB, Andrade JS, Herrmann HJ. 2006 Distribution of local fluxes in diluted porous media. Phys. Rev. E 74, 010401. (doi:10.1103/PhysRevE.74.010401) 21. Almeida MP, Maksean HA, Stanley HE, Andrade JS, Costa UMS. 1999 Inertial effects on fluid flow through disordered porous media. Phys. Rev. E 82, 5249. (doi:10.1103/PhysRevLett. 82.5249) 22. ANSYS FLUENT 12.0. 2009 Theory guide. Lebanon, NH: ANSYS Inc. 23. Munson BR, Young DF, Okiishi TH, Huebsch. WW 2009 Fundamentals of fluid mechanics. New York, NY: Wiley. 24. Saffman PG, Taylor GI. 1958 The penetration of a fluid into a porous medium or Hele– Shaw cell containing a more viscous liquid. Proc. R. Soc. Lond. A 245, 312–329. (doi:10.1098/ rspa.1958.0085)

Geometric packing constraints in egg phosphatidylcholine vesicles.

Geometric quantification of features in large flow fields.

Granular flow through an aperture: influence of the packing fraction.

Fluid flow in soft-walled tubes part 1: steady flow.

Piezoelectric energy harvesting in internal fluid flow.

Active polar fluid flow in finite droplets.

Air-fluid levels in the sphenoid sinus in epistaxis and nasal packing.

Cerebrospinal fluid flow. III. Pathological cerebrospinal fluid pulsations.

Flow visualization with blood analog fluid.

Microscale imaging of cilia-driven fluid flow.

Numerical flow studies in human carotid artery bifurcations: basic discussion of the geometric factor in atherogenesis.

Superconfinement tailors fluid flow at microscales.

Geometric pumping induced by shear flow in dilute liquid crystalline polymer solutions.

Effect of fluid friction on interstitial fluid flow coupled with blood flow through solid tumor microvascular network.

Mathematically modeling fluid flow and fluid shear stress in the canaliculi of a loaded osteon.

Basic cerebrospinal fluid flow patterns in ventricular catheters prototypes.

Flow cytometry in cerebrospinal fluid-rational use of laboratory services.

Accurate on-line mass flow measurements in supercritical fluid chromatography.

Numerical investigation of fluid flow in a chandler loop.

Cerebrospinal fluid flow dynamics in children with external ventricular drains.

Colonization, competition, and dispersal of pathogens in fluid flow networks.

Fluid Structural Analysis of Urine Flow in a Stented Ureter.

Peritumoral interstitial fluid flow velocity predicts survival in cervical carcinoma.

Fluid flow-induced calcium response in osteoclasts: signaling pathways.