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

0.98 MB·English

by  Alejandro David Mariotti

Papers in Physics, vol. 1, art. 010007 (2009)
www.papersinphysics.org
Received: 14 October 2009, Accepted: 28 December 2009
Edited by: A. C. Mart´ı
Licence: Creative Commons Attribution 3.0
DOI:10.4279/PIP.010007
ISSN 1852-4249
0
1 Parametric study of the interface behavior between two immiscible
0
2 liquids flowing through a porous medium

n
a
J Alejandro David Mariotti,1∗ Elena Brandaleze,2† Gustavo C. Buscaglia3‡

3
1
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,
]
n tilts and bends. The interface profiles depend on the physical properties of the liquids
y and on the velocity at which they are extracted. If the drainage flow rate, the liquids
d volumefraction in thedrainage flowandthephysicalproperties oftheliquidsareknown,
- the interface angle in the immediate vicinity of the outlet (θ) can be determined. In this
u
work, we define four nondimensional parameters that rule the fluid dynamical problem
l
f and, by means of a numerical parametric analysis, an equation to predict θ is developed.
.
s The equation is verified through several numerical assessments in which the parameters
c
are modified simultaneously and arbitrarily. In addition, the qualitative influence of each
i
s nondimensional parameter on theinterface shape is reported.
y
h
p
[ I. Introduction and thereby stimulate production. Normally, just

a small percentage of the oil in a reservoir can be
1
The fluid dynamics of the flow of two immiscible extracted, but water injection increases that per-
v
7 liquids through a porous medium plays a key role centage and maintains the production rate of the
9 in several engineering processes. Usually, though reservoir over a longer period of time. The water
1 the interest is focused on the extraction of one of displacestheoilfromthereservoirandpushesitto-
2 theliquids,thesimultaneousextractionofbothliq- wardsanoilproductionwell[1]. Inthe steelindus-
.
1 uids is necessary. This is the case of oilproduction try, this multiphase phenomenon occurs inside the
0 and of ironmaking. The water injection method blast furnace hearth, in which the porous medium
0
used in oil production consists of injecting water consists of coke particles. The slag and pig iron
1
backintothereservoir,usuallytoincreasepressure are stratified in the hearth and, periodically, they
:
v are drained through a lateral orifice. The under-
i ∗E-mail: [email protected]
X standingofthisflowiscrucialfortheproperdesign
†E-mail: [email protected]
r ‡E-mail: [email protected] andmanagementoftheblastfurnacehearth[2]. In
a bothexamplesabove,whentheliquidsaredrained,
1 Instituto Balseiro, 8400 San Carlos de Bariloche, Ar- a complex flow takes place in which the interface
gentina. between the liquids moves, tilts and bends.
2 Departamento de Metalurgia, Universidad Tecnolo´gica
Numerical simulation of multiphase flows in
Nacional Facultad Regional San Nicol´as, 2900 San
porous media is focused mainly in upscaling meth-
Nicol´as,Argentina.
ods, aimed at solving for large scale features of in-
3 Instituto de Ciˆencias Matema´ticas e de Computa¸ca˜o,
terest in such a way as to model the effect of the
UniversidadedeS˜aoPaulo,13560-970S˜aoCarlos,Brasil.
small scale features [3–5]. Other authors [6–8] use
010007-1
Papers in Physics, vol. 1, art. 010007 (2009) / A. D. Mariotti et al.
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 medium
when they are drained through an opening. The
effect of gravityon 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
predicttheinterfacetiltinthevicinityoftheorifice
(θ) 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.
Figure 1: Sketch of the numerical 2D domain.
II. Parametric Study
Thenumericalstudiesinthisworkwerecarriedout ratio, ρR = ρ1/ρ2, and nondimensional velocity,
by means of the program FLUENT 6.3.26. Differ- VR = V0ρ2L/µ2; where V0 is the outlet velocity
ent models to simulate the two-dimensional para- and L a reference length.
metricstudywereused. Thevolumeoffluid(VOF)
method was chosen to treat the interface problem i. Domain description
[9]. The drag force in the porous medium was
The numerical domain considered to carry out the
modeled by means ofthe sourceterm suggestedby
Forchheimer [10]. The source term for the ith mo- parametric study was a two-dimensional one com-
posed by the porous medium sub-domain and the
mentum equation is:
outlet sub-domain. The porous sub-domain is a
rectangle 10m wide and 10m tall. Inside of it, a
µ 1 −→ rigid, isotropic and homogeneous porous medium
S =− V + ρC|V |V . (1)
i (cid:18)α i 2 i(cid:19) was arranged. We use a porosity and particle di-
ameter of 0.32 and 0.006m, respectively.
FortheconstantsαandC inEq. (1),weusethe For the outlet domain we use a rectangle 0.02m
values proposed by Ergun [10]: 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
ε3d2
α= , (2) of the outlet domain is used to impose the outlet
150(1−ε)2
velocity. Figure 1 shows a complete description of
the domain.
A quadrilateral mesh with 2.2 × 104 cells was
(1−ε)
C =1.75 , (3) used, where the outlet sub-domain mesh consists
ε3d
of 200 elements in all the cases studied.
where ε is the porosity, d is the particle equivalent As boundary conditions, on edge 1 we define a
diameter, ρ is the density, V is the velocity and µ zerogaugepressureconditionnormaltothebound-
is the dynamic molecular viscosity. aryandimposethatonlythefluid1canentertothe
Consideringthat the subscript1 and2 represent domainthroughit. Onedge2the boundarycondi-
thefluid1andthefluid2respectively,threenondi- tionis the sameas onedge 1 butthe fluid consider
mensionalparameterswere consideredin the para- in this case is fluid 2. On edge 7 we impose a zero
metric study: viscosity ratio, µ = µ /µ , density gauge pressure normal to the boundary but in this
R 1 2
010007-2
Papers in Physics, vol. 1, art. 010007 (2009) / A. D. Mariotti et al.
Figure 2: Interface evolution when the initial posi- Figure 3: Interface evolution when the initial posi-
tion is below the outlet level, without gravity. tion is below the outlet, with gravity.
case the fluids can only leave the domain. On the
other edges (edges 4, 5, 6, and 8) we impose a wall
conditionwherethenormalandtangentialvelocity
is zero except for edge 3, at which the tangential
velocityisfreeandthestresstangentialtotheedge
is zero.
ii. Interface evolution
Toillustratehowaninterfacereachesthestationary
position froman initially horizontalone, 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-
Figure 4: Interface evolution when the initial posi-
facemodifiesitstilttoreachtheexitanditchanges
tion is above the outlet, with gravity.
its shape to reach the stationary profile.
Figures 3 and 4 show the interface evolution
whenthe gravityispresentbut the interfaceinitial
Two sets of curves were obtained, one consid-
positionis belowandabovethe outletlevelrespec-
ering the effect of gravity and the other without
tively.
considering it. Figure 5 shows the stationary in-
terface profiles for the different values of viscosity,
iii. Viscosity effect
without gravity. A value of V = 1.5× 104 and
R
Oneofthemostimportantparameterstomodifyis ρR = 1.17 were chosen. It is possible to observe
thedynamicalviscosityoffluid1. Wemaintainthe that, when fluid 1 has a viscosity higher than that
properties of the fluid 2 as the properties of water of fluid 2, the interface profile is above the outlet
(density 998Kg/m3, and dynamical viscosity 0.001 and points downwards at the outlet. If fluid 1 has
Pa.s)andthedensityoffluid1asthedensityofthe a lower viscosity, the opposite happens.
oil (850 Kg/m3). The dynamical viscosity of fluid When considering gravity, the value of V was
R
1wasvariedfromvaluessmallerthanthoseoffluid changed to 1×105 (V =10m/s), since for smaller
0
2, to values much greater. valuestheinterfacemaynotreachtheoutlet(thisis
010007-3
Papers in Physics, vol. 1, art. 010007 (2009) / A. D. Mariotti et al.
Figure 5: Stationary interface profiles modifying Figure 7: Stationary interface profiles for several
the fluid 1 viscosity without the gravity effect. values of VR, without gravity, for µR > 1 (µR =
35).
Figure 6: Stationary interface profiles modifying
the fluid 1 viscosity with the gravity effect. Figure 8: Stationary interface profiles for several
values of V , without gravity, for µ < 1 (µ =
R R R
0.01).
laterstudiedinFig. 10). Figure6showsthecurves
obtainedforthissituation,wheretheinterfaceonly
lies over the outlet level for the higher µ values.
R tained.
Figure7showstheeffectofV whentheviscosity
R
iv. Outlet velocity effect
of fluid 1 is greater than that of fluid 2, without
V0isvariedfromasmallvalue,similartotheporous gravity. It is possible to see that, as VR increases,
medium velocity (V0 =0.2m/s or VR =2000),to a the interface tilt at the outlet is maximal for VR =
verylargeone(V =50m/sorV =5×105). Main- 1×105.
0 R
taining the properties of fluid 2 similar to those of On the other hand, Fig. 8 shows the interface
water, two sets of curves were obtained (µ > 1 profileswhentheviscosityoffluid1issmallerthan
R
and µR <1), shown in Figs. 7 and 8, respectively. thatoffluid2. WeobservethatasVR increasesthe
When the effect of gravity was considered, two interface tends to the horizontal position.
additional sets of curves (Figs. 9 and 10) were ob- Figure 9 shows the different stationary interface
010007-4
Papers in Physics, vol. 1, art. 010007 (2009) / A. D. Mariotti et al.
Figure 9: Stationary interface profiles for several Figure 11: Stationary interface profiles for several
values of VR, with gravity,for µR >1 (µR =35). values ofthe density offluid 1, withVR =1.54 and
µ =35.
R
fluid 2 with the properties of water and the viscos-
ity of fluid 1 as 0.035Pa.s (µ = 35), the density
R
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 ρ increases, the interface profile
R
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 nondi-
Figure 10: Stationary interface profiles for several
mensional parameter on the interface behavior, an
values of V , with gravity,for µ <1 (µ =0.01).
R R R
equationthatpredictstheinterfaceangleattheim-
mediate vicinity of the outlet (θ) was crafted. For
positionswhenthegravityeffectispresentforµ > practicalreasons,thecaseswheregravityispresent
R
1. The effect of gravity is quite significant, the were considered to develop the equation.
interface ascends but only for the highest value of In Sect. 2.2, it is possible to see that when the
VR it lies above the outlet level. nondimensionalparametersρR,µR andVR arecon-
Figure 10showsthe curveswhenthe viscosityof stanttheinterfacechangesitsshapeuntilitreaches
fluid 1 is lower than that of fluid 2 (µR < 1). The a stationary profile.
behavior is different from that without gravity. In
For this reason,a fourthnondimensionalparam-
fact, there exists a minimum outlet velocity below
eter is considered, the volume fraction of fluid 1 in
which the interface does not reach the outlet.
the outlet flow (VF).
A generic expression [(Eq. (4)], consisting of
v. Density effect
three terms and containing 22 constants, was ad-
The effect ofthe density ratio onthe interface pro- justed by trial and error until satisfactory agree-
file wasstudied inthe presence ofgravity. Keeping ment with the numerical results was found.
010007-5
Papers in Physics, vol. 1, art. 010007 (2009) / A. D. Mariotti et al.
a -29.1 i 1.162×107 p 0.1 Case Fluid1 Fluid2 µ1 µ2 ρ1 ρ2
1 Heavy Water 0.4 0.005 850 998
b -0.04 j -0.18 q 0.72
oil solu-
c 0.1 k -1.64 r 2763.1 tion
2 Light Water 0.012 0.06 850 998
d 0.45 l -1 s 0.4
oil emul-
e 206 m 0.63 t -0.6 sion
3 – – 0.08 0.008 4000 4680
f 0.035 n 12.74 u -1.82 4 Kero- Water 24×10−4 0.001 780 998
g -0.05 o 0.085 v 3.2 sene
h -1.26 5 Ace- Water 3.3×10−4 0.001 791 998
tone
6 – – 0.003 0.01 400 720
7 – – 0.2 0.01 1500 2700
8 – – 0.3 0.004 3300 5940
Table 1: Constantvalues inthe genericexpression.
9 Light Hot 0.02 0.001 2800 7000
slag pig
iron
10 – – 0.5 0.013 500 1250
PT D ε 1/α C Resistance
A 0.005 0.3 10.88×107 1.81×104 Veryhigh 11 Mslaegdium piriogn 0.4 0.005 2800 7000
B 0.006 0.32 5.88×107 1.21×104 High 12-16 – – 0.08 0.008 4000 4680
C 0.02 0.2 3×107 1.75×104 Medium 17-21 Heavy pig 0.4 0.005 2800 7000
D 0.02 0.25 1.35×107 8400 Low slag iron
E 0.05 0.17 8.41×106 1.18×104 Verylow
Table 3: Cases description.
Table 2: Porous medium types.
Case PT ρR µR VR VF
1 B 1.17 80 5×104 2.8
2 B 1.17 0.2 5×104 87.2
θ =αµbRVRcρdR 3 B 1.17 10 10×104 28.4
+eµfVgρh exp(−iµj Vkρl VFm) 4 B 1.28 2.4 5×104 96.1
R R R R R R 5 B 1.26 0.33 1×105 96.8
+nµoRVRpρqRexp(−rµsRVRtρuRVFv) (4) 6 B 1.8 0.3 1×105 93.2
7 B 1.8 20 1×105 16.1
Table 1 shows the values of the constants in the 8 B 1.8 80 8×104 18.4
generic expression. 9 B 2.5 20 1×105 100.0
10 B 2.5 40 1×105 5.5
i. Equation verification 11 B 2.5 80 5×104 70.8
12 A 1.2 10.0 1×105 23.1
To verify that the generic expression (4) predicts 13 B 1.2 10.0 1×105 28.4
the value of θ correctly when the parameters are 14 C 1.2 10.0 1×105 41.5
arbitrarily modified, 21 additional numerical cases 15 D 1.2 10.0 1×105 53.4
were simulated. These cases cover a wide range of 16 E 1.2 10.0 1×105 63.7
physical properties of the liquids and of the char- 17 A 2.5 80.0 5×104 15.8
acteristics of the porous medium (by means of the 18 B 2.5 80.0 5×104 24.3
coefficients 1/α and C). 19 C 2.5 80.0 5×104 43.2
20 D 2.5 80.0 5×104 70.8
The interface angle obtained from the generic
21 E 2.5 80.0 5×104 94.0
expression (θ ) was compared with the interface
GE
angle obtained from the simulations (θ ). Table 2
S
showsthefiveporousmediumtypes(PT)thatwere
Table 4: Parameter values and PT for all cases
chosen.
described in Table 3.
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. parameterswascovered. Table3showsthedescrip-
For the remaining cases, the properties of the liq- tion of each numerical case, while Table 4 shows
uids were fixed at arbitrary values (fictitious liq- the corresponding nondimensional parameter val-
uids) so that a wide range of the nondimensional ues and PT.
010007-6
Papers in Physics, vol. 1, art. 010007 (2009) / A. D. Mariotti et al.
Acknowledgements - A. D. M. and E. B. are
gratefulforthesupportfromMetallurgicalDepart-
ment and DEYTEMA (UTNFSRN). G. C. B. ac-
knowledges partial financial support from CNPq
and FAPESP (Brazil).
[1] WCLyons,GJPlisga,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).
Figure 12: Comparison between the predictions of
the generic expressionand the numerical result for
[3] A Westhead, Upscaling for two-phase flows in
the 21 validation cases.
porous media, PhDthesis: CaliforniaInstitute
of Technology, Pasadena,California (2005).
We define an error (e = 100|∆θ|/180) as the [4] REEwing,TheMathematics ofreservoirsim-
percentage of the absolute value of the difference ulation, SIAM, Philadelphia (1983).
between the interface angles (∆θ = θ −θ ) di-
S GE
vided by the interface angle range (180◦). Figure [5] M A Cardoso, L J Durlofsky, Linearized
12 shows the comparison between the generic ex- reduced-order models for subsurface flow sim-
pression and the numerical cases. It is seen that ulation, J. Comput. Phys. 229, 681 (2010).
thegenericexpression(4)predictstheinterfacean-
[6] M J Blunt, Flow in porous media - pore-
gle, for the cases used in this study, with an error
network models and multiphase flow, Curr.
smaller than 10%.
Opin. Colloid Interface Sci. 6, 197 (2001).
[7] Z Chen,GHuan,YMa,Computational meth-
IV. Conclusions ods for multiphase flows in porous media,
SIAM, Philadelphia (2006).
A numerical study of the macroscopic interface
[8] Y Efendiev, T Houb, Multiscale finite element
behavior between two immiscible liquids flowing
methods for porous media flows and their ap-
through a porous medium, when they are drained
plications, Appl. Num. Math. 57, 577 (2007).
through an opening, has been reported. Four
nondimensional parameters that rule the fluid- [9] FLUENT6.3User’sGuide,FluentInc.(2006).
dynamical problem were identified. Thereby, a nu-
merical parametric analysis was developed where [10] J Bear, Dynamics of fluids in porous media,
the qualitative observation of the resulting inter- Dover Publications Inc., New York (1988).
faceprofilescontributestotheunderstandingofthe
effect of each parameter. In addition, a generic ex-
pressiontopredicttheinterfaceangleintheimme-
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
liquidsandporousmedia,thepredictionofθwithin
an error of 10% is considered satisfactory.
010007-7