Parametric study of the interface behavior between two immiscible liquids flowing through a porous medium

When two immiscible liquids that coexist inside a porous medium 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 (theta) 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 theta 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 liquids 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 percentage 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 towards an oil production well [1]. In the steel industry, 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 understanding 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 methods, aimed at solving for large scale features of interest in such a way as to model the effect of the small scale features [3][4][5]. Other authors [6][7][8] use the numerical methods to model the complex multiphase flow that takes place at the pore scale.
In this work, we numerically study the macroscopic behavior of the interface between two immiscible liquids flowing through a porous medium 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 orifice (θ) 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. Different models to simulate the two-dimensional parametric study were used. The volume of fluid (VOF) method was chosen to treat the interface problem [9]. The drag force in the porous medium was modeled by means of the source term suggested by Forchheimer [10]. The source term for the i th momentum equation is: For the constants α and C in Eq. (1), we use the values proposed by Ergun [10]: where ε is the porosity, d is the particle equivalent diameter, ρ is the density, V is the velocity and µ 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, µ R = µ 1 /µ 2 , density ratio, ρ R = ρ 1 /ρ 2 , and nondimensional velocity, where V 0 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 composed by the porous medium sub-domain and the outlet sub-domain. The porous sub-domain is a rectangle 10m wide and 10m tall. Inside of it, a rigid, isotropic and homogeneous porous medium 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 = 0.01m 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 × 10 4 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.

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 interface 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 respectively.

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 3 , and dynamical viscosity 0.001 Pa.s) and the density of fluid 1 as the density of the oil (850 Kg/m 3 ). The dynamical viscosity of fluid 1 was varied from values smaller than those of fluid 2, to values much greater.  Two sets of curves were obtained, one considering the effect of gravity and the other without considering it. Figure 5 shows the stationary interface profiles for the different values of viscosity, without gravity. A value of V R = 1.5 × 10 4 and ρ R = 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.
When considering gravity, the value of V R was changed to 1 × 10 5 (V 0 = 10m/s), since for smaller values the interface may not reach the outlet (this is  later 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 µ R values.
iv. Outlet velocity effect V 0 is varied from a small value, similar to the porous medium velocity (V 0 = 0.2m/s or V R = 2000), to a very large one (V 0 = 50m/s or V R = 5×105). Maintaining the properties of fluid 2 similar to those of water, two sets of curves were obtained (µ R > 1 and µ R < 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 ob-  tained. Figure 7 shows the effect of V R when the viscosity of fluid 1 is greater than that of fluid 2, without gravity. It is possible to see that, as V R increases, the interface tilt at the outlet is maximal for V R = 1 × 10 5 .
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 V R increases the interface tends to the horizontal position. Figure 9 shows the different stationary interface  positions when the gravity effect is present for µ R > 1. The effect of gravity is quite significant, the interface ascends but only for the highest value of V R it lies above the outlet level. Figure 10 shows the curves when the viscosity of fluid 1 is lower than that of fluid 2 (µ R < 1). The behavior is different from that without gravity. In fact, there exists a minimum outlet velocity below which the interface does not reach the outlet.

v. Density effect
The effect of the density ratio on the interface profile was studied in the presence of gravity. Keeping fluid 2 with the properties of water and the viscosity of fluid 1 as 0.035Pa.s (µ R = 35), the density of fluid 1 was varied from its original value to one three times smaller than that of fluid 2. In Fig. 11 it is seen that as ρ R 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 immediate vicinity of the outlet (θ) 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 ρ R , µ R and V R are constant 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 adjusted by trial and error until satisfactory agreement with the numerical results was found.   Table 1 shows the values of the constants in the generic expression.

i. Equation verification
To verify that the generic expression (4) predicts the value of θ correctly when the parameters are arbitrarily modified, 21 additional numerical cases were simulated. These cases cover a wide range of physical properties of the liquids and of the characteristics of the porous medium (by means of the coefficients 1/α and C).
The interface angle obtained from the generic expression (θ GE ) was compared with the interface angle obtained from the simulations (θ S ). Table 2 shows the five porous medium types (PT) that were chosen.
Some cases (1,2,4,5,9,11,(17)(18)(19)(20)(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 fixed at arbitrary values (fictitious liquids) so that a wide range of the nondimensional   Table 4: Parameter values and PT for all cases described in Table 3.
parameters was covered. Table 3 shows the description of each numerical case, while Table 4 shows the corresponding nondimensional parameter values and PT. We define an error (e = 100|∆θ|/180) as the percentage of the absolute value of the difference between the interface angles (∆θ = θ S − θ GE ) divided by the interface angle range (180 • ). Figure  12 shows the comparison between the generic expression 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 fluiddynamical problem were identified. Thereby, a numerical parametric analysis was developed where the qualitative observation of the resulting interface profiles contributes to the understanding of the effect of each parameter. In addition, a generic expression to predict the interface angle in the immediate vicinity of the outlet opening (θ) was developed. To verify that the generic equation predicts the value of θ correctly, 21 numerical cases with widely different parameters were simulated. Considering that the cases encompass a large class of liquids and porous media, the prediction of θ within an error of 10% is considered satisfactory.