Parametric study of the interface behavior between two immiscible liquids flowing through a porous médium

 

Alejandro David Mariotti,¹* Elena Brandaleze,² Gustavo C. Buscaglia³

* E-mail: mariotti.david@gmail.com
E-mail: ebrandaleze@frsn.utn.edu.ar
E-mail: gustavo.buscaglia@icmc.usp.br

Instituto Balseiro, 8400 San Carlos de Bariloche, Argentina.
Departamento de Metalurgia, Universidad Tecnológica Nacional Facultad Regional San Nicolás, 2900 San Nicolás, Argentina.
Instituto de Ciencias Matemáticas e de Computacáo, Universidade de Sao Paulo, 13560-970 Sao Carlos, Brasil.

 

When two immiscible liquids that coexist inside a porous médium are drained through an opening, a complex flow takes place in which the interface between the liquids moves, tilts and bends. The interface profiles depend on the physical properties of the liquids and on the velocity at which they are extracted. If the drainage flow rate, the liquids volume fraction in the drainage flow and the physical properties of the liquids are known, the interface angle in the immediate vicinity of the outlet (0) can be determined. In this work, we define four nondimensional parameters that rule the fluid dynamical problem and, by means of a numerical parametric analysis, an equation to predict 6 is developed. The equation is verified through several numerical assessments in which the parameters are modified simultaneously and arbitrarily. In addition, the qualitative influence of each nondimensional parameter on the interface shape is reported.

I. Introduction

The fluid dynamics of the flow of two immiscible liquids through a porous medium plays a key role in several engineering processes. Usually, though the interest is focused on the extraction of one of the liquids, the simultaneous extraction of both liq-uids is necessary. This is the case of oil production and of ironmaking. The water injection method used in oil production consists of injecting water back into the reservoir, usually to increase pressure and thereby stimulate production. Normally, just a small percentage of the oil in a reservoir can be extracted, but water injection increases that per-centage and maintains the production rate of the reservoir over a longer period of time. The water displaces the oil from the reservoir and pushes it to-wards an oil production well [1]. In the steel indus-try, this multiphase phenomenon occurs inside the blast furnace hearth, in which the porous medium consists of coke particles. The slag and pig iron are stratified in the hearth and, periodically, they are drained through a lateral orifice. The under-standing of this flow is crucial for the proper design and management of the blast furnace hearth [2]. In both examples above, when the liquids are drained, a complex flow takes place in which the interface between the liquids moves, tilts and bends.

Numerical simulation of multiphase flows in porous media is focused mainly in upscaling meth-ods, aimed at solving for large scale features of interest in such a way as to model the effect of the small scale features [3-5]. Other authors [6-8] use the numerical methods to model the complex mul-tiphase flow that takes place at the pore scale.

In this work, we numerically study the macro-scopic behavior of the interface between two im-miscible liquids flowing through a porous médium when they are drained through an opening. The effect of gravity on this phenomenon is considered. We define four nondimensional parameters that rule the fluid dynamical problem and, by means of a numerical parametric analysis, an equation to predict the interface tilt in the vicinity of the orífice (9) is developed. The equation is verified through several numerical cases where the parameters are varied simultaneously and arbitrarily. In addition, the qualitative influence of each non-dimensional parameter on the interface shape is reported.

II. Parametric Study

The numerical studies in this work were carried out by means of the program FLUENT 6.3.26. Differ-ent models to simúlate the two-dimensional parametric study were used. The volume of fluid (VOF) method was chosen to treat the interface problem [9]. The drag forcé in the porous médium was modeled by means of the source term suggested by Forchheimer [10]. The source term for the i mo-mentum equation is:

Figure 1: Sketch of the numerical 2D domain.

where e is the porosity, d is the particle equivalent diameter, p is the density, V is the velocity and p, is the dynamic molecular viscosity.

Considering that the subscript 1 and 2 represent the fluid 1 and the fluid 2 respectively, three nondimensional parameters were considered in the parametric study: viscosity ratio, p,R = p,\/p,2, density ratio, pR = pi/p2, and nondimensional velocity Vr = V0P2L/P2', where Vq is the outlet velocity and L a reference length.

i. Domain description

The numerical domain considered to carry out the parametric study was a two-dimensional one com-posed by the porous médium sub-domain and the outlet sub-domain. The porous sub-domain is a rectangle lOm wide and lOm tall. Inside of it, a rigid, isotropic and homogeneous porous médium was arranged. We use a porosity and particle diameter of 0.32 and 0.006m, respectively.

For the outlet domain we use a rectangle 0.02m wide and with a height L = O.Olm divided into two equal parts and located at the center of one of the lateral edges. The part located at the end of the outlet domain is used to impose the outlet velocity. Figure 1 shows a complete description of the domain.

A quadrilateral mesh with 2.2 x 10⁴ cells was used, where the outlet sub-domain mesh consists of 200 elements in all the cases studied.

As boundary conditions, on edge 1 we define a zero gauge pressure condition normal to the boundary and impose that only the fluid 1 can enter to the domain through it. On edge 2 the boundary condition is the same as on edge 1 but the fluid consider in this case is fluid 2. On edge 7 we impose a zero gauge pressure normal to the boundary but in this case the fluids can only leave the domain. On the other edges (edges 4, 5, 6, and 8) we impose a wall condition where the normal and tangential velocity is zero except for edge 3, at which the tangential velocity is free and the stress tangential to the edge is zero.

 

Figure 2: Interface evolution when the initial posi-tion is below the outlet level, without gravity.

 

ii. Interface evolution

To illustrate how an interface reaches the stationary position from an initially horizontal one, three sets of curves were obtained.

Figure 2 shows the interface evolution for the case without the gravity effect and the interface initial position is below the outlet level. The inter-face modifies its tilt to reach the exit and it changes its shape to reach the stationary profile.

Figures 3 and 4 show the interface evolution when the gravity is present but the interface initial position is below and above the outlet level respec-tively.

iii. Viscosity effect

One of the most important parameters to modify is the dynamical viscosity of fluid 1. We maintain the properties of the fluid 2 as the properties of water (density 998Kg/m³, and dynamical viscosity 0.001 Pa.s) and the density of fluid 1 as the density of the oil (850 Kg/m³). The dynamical viscosity of fluid 1 was varied from values smaller than those of fluid 2, to values much greater.

 

Figure 3: Interface evolution when the initial posi-tion is below the outlet, with gravity.

Figure 4: Interface evolution when the initial posi-tion is above the outlet, with gravity.

Two sets of curves were obtained, one considering the effect of gravity and the other without considering it. Figure 5 shows the stationary in-terface profiles for the different valúes of viscosity, without gravity. A valué of Vr = 1.5 x 10⁴ and Pr = 1.17 were chosen. It is possible to observe that, when fluid 1 has a viscosity higher than that of fluid 2, the interface profile is above the outlet and points downwards at the outlet. If fluid 1 has a lower viscosity the opposite happens.

 

Figure 5: Stationary interface profiles modifying the fluid 1 viscosity without the gravity effect.

Figure 6: Stationary interface profiles modifying the fluid 1 viscosity with the gravity effect.

When considering gravity the valué of Vr was changed to 1 x 10⁵ (Vb = lOm/s), since for smaller valúes the interface may not reach the outlet (this islater studied in Fig. 10). Figure 6 shows the curves obtained for this situation, where the interface only lies over the outlet level for the higher ¡ir valúes.

iv. Outlet velocity effect

Vo is varied from a small valué, similar to the porous médium velocity (Vo = 0.2m/s or Vr = 2000), to a very large one (Vo = 50m/s or Vr = 5x 105). Main-taining the properties of fluid 2 similar to those of water, two sets of curves were obtained (p,R > 1 and ¡ir < 1), shown in Figs. 7 and 8, respectively. When the effect of gravity was considered, two additional sets of curves (Figs. 9 and 10) were obtained.

Figure 7: Stationary interface profiles for several valúes of Vr, without gravity, for ¡ir > 1 (p,R = 35).

Figure 8: Stationary interface profiles for several valúes of Vr, without gravity, for ¡ir < 1 (p,R = 0.01).

 

Figure 7 shows the effect of Vr when the viscosity of fluid 1 is greater than that of fluid 2, without gravity. It is possible to see that, as Vr increases, the interface tilt at the outlet is maximal for Vr = 1 x 10⁵.

On the other hand, Fig. 8 shows the interface profiles when the viscosity of fluid 1 is smaller than that of fluid 2. We observe that as Vr increases the interface tends to the horizontal position.

Figure 9 shows the different stationary interface positions when the gravity effect is present for ¡ir > 1. The effect of gravity is quite significant, the interface ascends but only for the highest valué of Vr it lies above the outlet level.

 

Figure 9: Stationary interface profiles for several valúes of Vr, with gravity, for ¡ir > 1 (jj,r = 35).

Figure 10: Stationary interface profiles for several valúes of Vr, with gravity for ¡ir < 1 (p,R = 0.01).

 

Figure 10 shows the curves when the viscosity of fluid 1 is lower than that of fluid 2 (p,R < 1). The behavior is different from that without gravity. In fact, there exists a mínimum outlet velocity below which the interface does not reach the outlet.

 

Figure 11: Stationary interface profiles for several valúes of the density of fluid 1, with Vr = 1.5⁴ and Hr = 35.

fluid 2 with the properties of water and the viscosity of fluid 1 as 0.035Pa.s (jj,r = 35), the density of fluid 1 was varied from its original valué to one three times smaller than that of fluid 2. In Fig. 11 it is seen that as pR increases, the interface profile ascends significantly with a less significant change in the tilt angle at the outlet.

III. Generic expression

From the study on the influence of each nondimensional parameter on the interface behavior, an equation that predicts the interface angle at the im-mediate vicinity of the outlet (9) was crafted. For practical reasons, the cases where gravity is present were considered to develop the equation.

In Sect. 2.2, it is possible to see that when the nondimensional parameters pR, p,R and Vr are con-stant the interface changes its shape until it reaches a stationary profile.

For this reason, a fourth nondimensional parameter is considered, the volume fraction of fluid 1 in the outlet flow (VF).

A generic expression [(Eq. (4)], consisting of three terms and containing 22 constants, was ad-justed by trial and error until satisfactory agree-ment with the numerical results was found.

Table 1: Constant values in the generic expression.

a	-29.1	i	1.162 x 10	p	0.1
b	-0.04	j	-0.18	q	0.72
c	0.1	k	-1.64	r	2763.1
d	0.45	1	-1	s	0.4
e	206	m	0.63	t	-0.6
f	0.035	n	12.74	u	-1.82
g	-0.05	o	0.085	V	3.2
h	-1.26

Table 2: Porous médium types.

PT	D	e	1/a	C	Resistance
A	0.005	0.3	10.88 x 10	1.81 x 10	Very high
B	0.006	0.32	5.88 x 10	1.21 x 10	High
C	0.02	0.2	3 x 10	1.75 x 10	Médium
D	0.02	0.25	1.35 x 10	8400	Low
E	0.05	0.17	8.41 x 10	1.18 x 10	Very low

 

9 = ap_RV_Rp_R

+ e^RVRPRe^XP(-^Íl^J'RVRP^lRVF^m)

+ np_RV_Rp^q_R exp(- rp_RV_Rp_RVF^v) (4)

Table 1 shows the valúes of the constants in the generic expression.

i. Equation verification

To verify that the generic expression (4) predicts the valué of 9 correctly when the parameters are arbitrarily modiñed, 21 additional numerical cases were simulated. These cases cover a wide range of physical properties of the liquids and of the char-acteristics of the porous médium (by means of the coefñcients 1/a and C).

The interface angle obtained from the generic expression (9ge) was compared with the interface angle obtained from the simulations (9s). Table 2 shows the ñve porous médium types (PT) that were chosen.

Some cases (1, 2, 4, 5, 9, 11, 17-21) were chosen based on the possible combinations of immiscible liquids that can be manipulated in real situations. For the remaining cases, the properties of the liquids were ñxed at arbitrary valúes (ñctitious liquids) so that a wide range of the nondimensional parameters was covered. Table 3 shows the description of each numerical case, while Table 4 shows the corresponding nondimensional parameter val-ues and PT.

Table 3: Cases description.

Case	Fluid 1	Fluid 2	Mi	A¿2	Pi	P2
1	Heavy oil	Water solu-tion	0.4	0.005	850	998
2	Light oil	Water emul-sion	0.012	0.06	850	998
3	-	-	0.08	0.008	4000	4680
4	Kerosene	Water	24 x 10-4	0.001	780	998
5	Ace-tone	Water	3.3 x 10-	0.001	791	998
6	-	-	0.003	0.01	400	720
7	-	-	0.2	0.01	1500	2700
8	-	-	0.3	0.004	3300	5940
9	Light slag	Hot pig iron	0.02	0.001	2800	7000
10	-	-	0.5	0.013	500	1250
11	Médium slag	pig iron	0.4	0.005	2800	7000
12-16	-	-	0.08	0.008	4000	4680
17-21	Heavy slag	pig iron	0.4	0.005	2800	7000

Table 4: Parameter values and PT for all cases described in Table 3.

Case	PT	fs R	/-*. R	V ^-¿		VF
i	B	1.17	80	5 x	10	2.8
2	B	1.17	0.2	5 x	10	87.2
3	B	1.17	10	10 x	10	28.4
4	B	1.28	2.4	5 x	10	96.1
5	B	1.26	0.33	1 X	10	96.8
6	B	1.8	0.3	1 X	10	93.2
7	B	1.8	20	1 X	10	16.1
8	B	1.8	80	8 x	10	18.4
9	B	2.5	20	1 X	10	100.0
10	B	2.5	40	1 X	10	5.5
11	B	2.5	80	5 x	10	70.8
12	A	1.2	10.0	1 X	10	23.1
13	B	1.2	10.0	1 X	10	28.4
14	C	1.2	10.0	1 X	10	41.5
15	D	1.2	10.0	1 X	10	53.4
16	E	1.2	10.0	1 X	10	63.7
17	A	2.5	80.0	5 x	10	15.8
18	B	2.5	80.0	5 x	10	24.3
19	C	2.5	80.0	5 x	10	43.2
20	D	2.5	80.0	5 x	10	70.8
21	E	2.5	80.0	5 x	10	94.0

 

Figure 12: Comparison between the predictions of the generic expression and the numerical result for the 21 validation cases.

 

We define an error (e = 1OO|A0|/18O) as the percentage of the absolute valué of the difference between the interface angles (A9 = 9,s - Oge) di-vided by the interface angle range (180°). Figure 12 shows the comparison between the generic ex-pression and the numerical cases. It is seen that the generic expression (4) predicts the interface angle, for the cases used in this study, with an error smaller than 10%.

IV. Conclusions

A numerical study of the macroscopic interface behavior between two immiscible liquids flowing through a porous medium, when they are drained through an opening, has been reported. Four nondimensional parameters that rule the fluid-dynamical problem were identified. Thereby, a nu-merical parametric analysis was developed where the qualitative observation of the resulting inter-face profiles contributes to the understanding of the effect of each parameter. In addition, a generic ex-pression to predict the interface angle in the imme-diate vicinity of the outlet opening (θ) was devel-oped. To verify that the generic equation predicts the value of θ correctly, 21 numerical cases with widely different parameters were simulated. Con-sidering that the cases encompass a large class of liquids and porous media, the prediction of θ within an error of 10% is considered satisfactory.

Acknowledgements - A. D. M. and E. B. are grateful for the support from Metallurgical Depart-ment and DEYTEMA (UTNFSRN). G. C. B. ac-knowledges partial financial support from CNPq and FAPESP (Brazil).

[1] W C Lyons, G J Plisga, Standard Handbook of Petroleum & Natural Gas Engineering - 2nd ed, Gulf Professional Publishing, Burlington (2005).

[2] A K Biswas, Principles of blast furnace iron-making: theory and practice, Cootha Publish-ing House, Brisbane (1981).

[3] A Westhead, Upscaling for two-phase flows in porous media, PhD thesis: California Institute of Technology, Pasadena, California (2005).

[4] R E Ewing, The Mathematics of reservoir sim-ulation, SIAM, Philadelphia (1983).

[5] M A Cardoso, L J Durlofsky, Linearized reduced-order models for subsurface flow sim-ulation, J. Comput. Phys. 229, 681 (2010).

[6] M J Blunt, Flow in porous media - pore-network models and multiphase flow, Curr. Opin. Colloid Interface Sci. 6, 197 (2001).

[7] Z Chen, G Huan, Y Ma, Computational meth-ods for multiphase flows in porous media, SIAM, Philadelphia (2006).

[8] Y Efendiev, T Houb, Multiscale finite element methods for porous media flows and their ap-plications, Appl. Num. Math. 57, 577 (2007).

[9] FLUENT 6.3 User’s Guide, Fluent Inc. (2006).

[10] J Bear, Dynamics of fluids in porous media, Dover Publications Inc., New York (1988).