Contents lists available at ScienceDirect

ADVANCING RESEARCH EVOLVING SCIENCE

Petroleum

journal homepage: www.keaipublishing.com/en/journals/petlm

Quantitative prediction of residual wetting film generated in mobilizing a two-phase liquid in a capillary model

CrossMark

Harsh Joshi, Liming Dai*

Industrial Systems Engineering, University of Regina, Regina, SK S4S 0A2, Canada

ARTICLE INFO

Article history: Received 24 August 2015 Received in revised form 15 October 2015 Accepted 16 October 2015

Keywords: Capillary model Residual wetting film Higher viscosity liquid Water flooding Numerical simulation Immiscible flow Oil slug mobilization

ABSTRACT

This research studies the motion of immiscible two-phase liquid flow in a capillary tube through a numerical approach employing the volume of fluid method, for simulating the core-annular flow and water flooding in oil reservoirs of porous media. More specifically, the simulations are a representation of water flooding at a pore scale. A capillary tube model is established with ANSYS Fluent and verified. The numerical results matches well with the existing data available in the literature. Penetration of a less viscous liquid in a liquid of higher viscosity and the development of a residual wetting film of the higher viscosity liquid are thoroughly investigated. The effects of Capillary number, Reynolds Number and Viscosity ratio on the residual wetting film are studied in detail, as the thickness is directly related to the residual oil left in the porous media after water flooding. It should be noticed that the liquids considered in this research can be any liquids of different viscosity not necessarily oil and water. The results of this study can be used as guidance in the field of water flooding.

Copyright © 2015, Southwest Petroleum University. Production and hosting by Elsevier B.V. on behalf of KeAi Communications Co., Ltd. This is an open access article under the CC BY-NC-ND

license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction

In this paper, we consider a very important immiscible multiphase flow case of a non-wetting liquid displacing a wetting liquid which subsequently leaves behind a thin film of wetting liquid. This kind of flow commonly occurs in a petroleum reservoirs and it is known as drainage in Petroleum Engineering. Another name for this type of flow is core-annular flow. This flow is mimicking what happens in an ideal oil reservoir. A pore of the reservoir maybe filled with oil and then water flooding is commenced to recover the oil from the reservoir. The water will enter the pore and displace some amount of oil while leaving behind some oil in a form of a thin film. Since a single pore is

* Corresponding author. E-mail address: liming.dai@uregina.ca (L Dai). Peer review under responsibility of Southwest Petroleum University.

being studied rather than the entire reservoir, this study can be considered a pore scale study.

A similar case of non-wetting gas phase displacing a wetting liquid has been widely studied and number of different equations and correlations has been proposed to calculate the wetting film thickness left behind. One of the first equation to predict the wetting film thickness was proposed as 5 = 0.5Ca0 5 by Ref. [1] and valid for capillary numbers from 10~4 to 10~2. Ref. [2] proposed an equation of the form 5 = 1.3Ca3 derived theoretically and tested the equation by conducting experiments for a range of capillary number. The equation is valid for the range of capillary numbers: 10~6 to 10~2. Gas—liquid experiments conducted by Ref. [3] predicted the wetting film thickness to be dependent on the gas bubble length and contradicted the conclusion of Ref. [2]. Ref. [4] has offered an explanation for discrepancy between experimental data and theoretical predictions based on the presence of trace amount of surfactants giving rise to Marangoni surface flows in the fluids used. The author still could not fully explain the experimental and theoretical discrepancy as noted by Ref. [5]. Gas displacing wetting liquid experiments has also been carried out by Refs. [6,7]. It is noted by the author of [2] that for values of capillary numbers larger than 3 x 10~3, the results of [1,2,6,7] are in broad

http://dx.doi.org/10.1016/j.petlm.2015.10.005

2405-6561/Copyright © 2015, Southwest Petroleum University. Production and hosting by Elsevier B.V. on behalf of KeAi Communications Co., Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

agreement with the film thickness equation proposed by Ref. [1]. Below the capillary number of 3 x 10~3, there is a wide spread disagreement and different equations to predict the wetting film thickness are proposed. In experiments to measure wetting film thickness conducted by Ref. [8], the film thickness values are over predicted because of the measurement techniques employed by the author. Authors of Ref. [9] have used computational Fluid Dynamics to calculate the thickness of wetting film in a gas displacing liquid flow. The simulation results of this study were based on the same assumption made by Ref. [2]. Therefore, it cannot confirm the results of Ref. [2] like an actual experiment. Experiments of [10] in vertical tubes also found discrepancies between experimental and theoretical values obtained from equations of Ref. [1,2] respectively. The theoretical values under predict the film thickness value. The author of Ref. [11] has used Finite element method to confirm the prediction of Ref. [9] that finite Reynolds number (inertial) effects have minor effect on wetting film thickness. The author used much higher Reynolds numbers than previous author.

There are no equations or correlations out there for the case of a non-wetting liquid displacing a wetting liquid. The experimental data is also very scarce and controversial. The equation proposed by Ref. [2] is sometimes used to calculate the film thickness for liquid—liquid flow. However, this is incorrect because the original equation derived by Ref. [2] was obtained with the assumption of gas displacing liquid. The authors of Refs. [12—14] have conducted experiments to measure the wetting film thickness in liquid—liquid flow which can be used to validate a numerical simulation results. Ref. [15] has proposed an asymptotic theory by using boundary-integral method for the wetting film but does not offer any testable prediction for film thickness. The latest experiment on the wetting film thickness in immiscible liquid—liquid flow is done by Ref. [5]. The author has tabulated his experimental data in form of a table and also included in the table are the experimental values of Refs. [12—14]. A testable wetting film thickness equation has been derived from his theory which can be used to fit the experimental data. They have predicted lower capillary number film thickness behavior to be proportional to Ca2 and higher capillary number film thickness to approach a constant value of 0.2. This study also improves on the previous simulations of the group members [16,17]. They have measured the residual wetting liquid thickness in a liquid—liquid flow but their comparison is done to equations for liquid thickness in gas—liquid flow. This approach does not accurately verify their results and this study aims to improve upon that. The lack of theoretical equations for

liquid—liquid flow leads to the necessity of the numerical simulations.

The purpose of this study is to numerically study the wetting film thickness in the liquid—liquid immiscible two phase flow and compare it to the scarce literature data. Numerical model makes it possible to visualize and thoroughly study the residual wetting film phenomena along with the residual wetting film thickness. Effect of viscosity ratio and Reynolds number on the residual wetting film thickness are also studied. To the best of our knowledge, we have not encountered any papers studying water flooding with VOF method. We have also not encountered any papers using VOF method to simulate immiscible water—oil core annular flow in a single pore and measuring the thickness of oil film left behind. The study is aimed at helping to understand and visualize the process of water flooding in a single pore as well as better understand mechanisms of oil recovery in the porous media.

2. Method

2.1. Mathematical modeling

The software employed to perform the simulations is a popular software by ANSYS known as Fluent. The multiphase simulation in Fluent is based on a very popular approach known as the Volume of Fluid method (VOF). VOF method is used to track the interface between two immiscible liquids so that a wetting film thickness can be measured and compared to published results.

Some of the major assumptions used for the simulations in this paper are: laminar flow, incompressible fluid, and gravity effects are negligible. The gravity effects are negligible because of the small diameter tube employed in the simulation which leads to a small bond number [5]. The governing equations of this flow are the conservation of mass and momentum, along with the volume fraction equation [18].

The cylindrical capillary tube is modeled by two-dimensional axisymmetric option. The purpose of using axisymmetric option is that it enables us to model only half of the actual flow domain, thus resulting in fewer elements and saving computational time. The schematic of the geometry along with the boundary conditions are shown in detail in the Fig. 1. As seen in the figure, the left end of the tube has a velocity-inlet condition where a fixed velocity is applied. The right end has a pressure-outlet condition is applied and a gauge pressure of 0 Pa is set meaning atmospheric pressure. At the wall, there is a no slip condition applied.

Fig. 1. The geometry and dimensions of the flow of being studied in this paper also known as core-annular flow.

Fig. 1 shows that the diameter of the cylindrical tube is 0.5 mm and the length is 5 mm. This tube is employed for all the simulations in this paper.

2.2. Solution methodology

The simulation is conducted by using commercial CFD software ANSYS Fluent to simulate two-dimensional, axisymmetric, transient, laminar, immiscible multiphase flow in a small diameter cylindrical capillary tube.

2.3. VOF method

The volume-of-fluid method is used to track the interface between two immiscible liquids by solving the volume fraction equations given above [19]. Courant number which is important for stability of the fluid flow is set to 0.25. The time-step in the simulation is then determined by using the variable time stepping method determined based on the set courant number, the mesh density and the cell velocity. The equation for Courant number is as follows

Co = VfluidDx (1)

where Ax is the element size, At is the time step, and Vflujd is the fluid velocity in the element. The VOF courant number for the VOF equation and the global courant number for rest of the transport equations are both set to 0.25.

2.4. Solver options

The multiphase volume-of-Fluid model in Fluent 13 is selected for the simulation of the multiphase flow. Implicit body force option is turned on to stabilize the solution as recommended in Fluent 13 user's guide [20]. A first order, Non-iterative time advancement method is selected for the transient formulation because it saves computational time according to Fluent 13.0 user's guide [20]. The other solution methods are chosen according to Ref. [21]. The pressure—velocity coupling is obtained by the PISO method. Green-Gauss node based is chosen for the gradient. PRESTO! is used for pressure discretization. For momentum equation, QUICK is chosen and for volume fraction equation, explicit Geo-Reconstruct is chosen. The explicit Geo-Reconstruct scheme works by using a piecewise-linear

approach to represent the interface [22]. The pressure and Momentum correction tolerances for Non-iterative solver controls are set to 0.01 and 0.001 respectively. It was ensured that the absolute values of residuals achieved were sufficiently low (O(10-7) for x and y velocities and O(10 10) for continuity). The contact angle for all the cases is set to 180° for the non-wetting fluid meaning the wetting fluid perfectly wets the tube wall. Fluent employs the continuum surface force (CSF) model proposed by Ref. [23] to model surface tension effects.

3. Results and discussion

A small diameter cylindrical tube is used to simulate multiphase immiscible liquid—liquid flow. The simulations are done to match the wetting film thickness obtained from using volume-of-fluid method and comparing them to the latest experimental published values by Ref. [5]. A small circular tube as mentioned earlier has been used to simulate the immiscible multiphase flow. Impacts of capillary number, Reynolds number, and viscosity ratio on the residual wetting thickness are studied to see how they impact the thickness of the residual wetting film. Simulations are very useful because of its capacity to simulate a flow with different parameters that might otherwise be hard to control in an experiment. It is difficult or impossible to determine relationships between wetting film thickness and viscosity ratio and Reynolds number via experiments and this is where numerical simulations come in handy.

3.1. Meshing

The meshing of the geometry was done in Ansys workbench. The quality of the mesh was judged by three different metrics in the Ansys Workbench: element quality, skewness, and orthogonal quality. The meshing used in the simulations is shown in Fig. 2. The mesh has 136,000 quadrilateral elements and 0.01 mm of length near the wall is meshed finer and the rest 0.24 mm length away from the wall is meshed coarse. The advantage of this type of meshing is that the thin films can be captured accurately without increasing the meshing elements in the rest of the domain. The quadrilateral elements in the fine region has a width of 1 mm and a length of 2.9412 mm. The quadrilateral elements in the coarse region has a width of 4.8 mm and a length of 2.9412 mm. The mesh satisfies the five element across the film rule [21]. A grid independence test is done to pick

Fig. 2. The mesh is fine near the wall to capture the thin films accurately.

out the mesh that can produce accurate results. See the author's thesis [26] for detailed description of the meshing used and the verification of the meshing.

mannular M-core

3.2. Immiscible fluid interface

It is a known fact that the interface between two immiscible liquids is very thin and on the order of an angstrom [24]. Also, no gradual change exists from one liquid to the other. The volume fraction changes of the liquids are sudden at the interface. The VOF method cannot capture the sharp interface and hence, the resulting interface will appear smeared rather than sharp. The smeared interface is very thin but not as thin as it would be in reality. This is why report an upper and lower bound for the wetting film thickness is reported in Table 1. The wetting film thickness can be anywhere between the upper and lower bound because of the smearing effect. The smearing effect is shown in Fig. 3. The lower and upper bound can be averaged and reported as one value and this approach is used in Tables 2 and 3.

3.3. Capillary number and film thickness

The three main dimensionless parameters used in this study are Capillary number, Reynolds number, and Viscosity ratio as defined below, namely

Ca — Vinlet(moil/Toil-water)

V is the velocity at the inlet of the tube, m is the viscosity of the oil, and g is the water—oil interfacial tension.

Poil $ Vinlet $D Moil

p is the density of the oil, V is the velocity at the inlet of the tube, D is the diameter of the tube, and m is the oil viscosity.

mannular is the viscosity of the outside annual liquid (oil) in the tube and mcore is the viscosity of the core liquid (water) in the middle surrounded by the annular liquid.

The wetting film thickness is proportional to the tube radius as observed from the Bretherton's equation in Ref. [2]. In order to match experimental results which may employ tubes of different diameters and simulations, a normalized film thickness is used. The normalized film is dimensionless and is obtained by dividing the actual film thickness, obtained from the simulation, by the radius of the tube employed in the simulation. The normalized film thickness will be denoted by the symbol 5. Despite Breh-terton's equation not being valid for the entire range of capillary numbers, as proposed by the author, it is generally believed to be accurate in the capillary range of 10~5 to 10~2 [5]. The equation is stated below, once again, for convenience:

Bretherton's Equation : 5 — 1.3Ca2/3

In the first case, a wetting film thickness, at a given capillary number is of interest. The capillary number and Reynolds number are calculated based upon the properties of the wetting liquid (oil) by using the equations (2) and (3) as shown in Table 1. A small inlet velocity can create a long computation time with which to obtain the solution because the flow velocity is very slow. Large inlet velocity can lead to a high courant number, resulting in a smaller time step and longer computation time. Hence, velocity range is kept between 0.000788 m/s and 0.082095 m/s, as shown in Table 1. Reynolds number is kept small (<2 for all cases) to avoid any inertial effect affecting the film thickness as shown in Table 1. The viscosity ratio is kept constant at 147.2 in all cases to ensure that it does not affect the wetting film thickness. Flow is simulated at various capillary numbers and constant viscosity ratios to verify its effect upon the residual wetting film thickness. The range of capillary numbers simulated is from 0.00096 to 2. The range is chosen because most microflow conditions

Table 1

Comparison of the simulations film thickness with the published experimental values. A dash is used to represent no value in that paper at that capillary number.

Ca p water Poil m water Moil Yoil-water M Vinlet Re LB UB [12] [13] [14] [5] Equation

(kg/ m3) (kg/ m3) (kg/(mS)) (kg/(mS)) (N/m) (m/s) (5) (5) (5) (5) (5) (5) (5) (5)

0.00096 998.2 860 0.001003 0.148 0.1215 147.6 0.0007881081 0.0022897736 0.015 0.018 — — — 0.023 0.013

0.01 998.2 860 0.001003 0.148 0.1215 147.6 0.0082094595 0.0238518079 0.044 0.076 0.059 — — 0.061 0.062

0.022 998.2 860 0.001003 0.148 0.1215 147.6 0.0180608108 0.0524739774 0.078 0.112 0.067 0.088 — — —

0.033 998.2 860 0.001003 0.148 0.1215 147.6 0.0270912162 0.078710966 0.100 0.131 — — 0.093 — —

0.048 998.2 860 0.001003 0.148 0.1215 147.6 0.0394054054 0.1144886779 0.121 0.158 — 0.130 — — —

0.05 998.2 860 0.001003 0.148 0.1215 147.6 0.0410472973 0.1192590394 0.123 0.164 — 0.120 — — —

0.072 998.2 860 0.001003 0.148 0.1215 147.6 0.0591081081 0.1717330168 0.151 0.188 — — — 0.130 —

0.075 998.2 860 0.001003 0.148 0.1215 147.6 0.0615709459 0.1788885592 0.158 0.189 — 0.150 — — —

0.082 998.2 860 0.001003 0.148 0.1215 147.6 0.0673175676 0.1955848247 0.160 0.200 — 0.170 0.140 — —

0.1 998.2 860 0.001003 0.148 0.1215 147.6 0.0820945946 0.2385180789 0.179 0.212 — 0.170 — — —

0.17 998.2 860 0.001003 0.148 0.003375 147.6 0.0038766892 0.0112633537 0.220 0.262 — 0.210 — — —

0.19 998.2 860 0.001003 0.148 0.003375 147.6 0.0043327703 0.0125884542 0.234 0.266 — 0.240 — — —

0.2 998.2 860 0.001003 0.148 0.003375 147.6 0.0045608108 0.0132510044 0.236 0.269 — 0.240 — 0.210 —

0.23 998.2 860 0.001003 0.148 0.003375 147.6 0.0052449324 0.015238655 0.245 0.284 — 0.230 — — —

0.25 998.2 860 0.001003 0.148 0.003375 147.6 0.0057010135 0.0165637555 0.254 0.286 — — — 0.210 —

0.26 998.2 860 0.001003 0.148 0.003375 147.6 0.0059290541 0.0172263057 0.255 0.289 — 0.300 — — —

0.29 998.2 860 0.001003 0.148 0.003375 147.6 0.0066131757 0.0192139564 0.260 0.302 — 0.240 — — —

0.5 998.2 860 0.001003 0.148 0.003375 147.6 0.011402027 0.033127511 0.295 0.336 — 0.270 — — —

0.88 998.2 860 0.001003 0.148 0.003375 147.6 0.0200675676 0.0583044193 0.323 0.360 — — — 0.270 —

1 998.2 860 0.001003 0.148 0.003375 147.6 0.0228040541 0.0662550219 0.331 0.363 — 0.320 — 0.270 —

2 998.2 860 0.001003 0.148 0.003375 147.6 0.0456081081 0.1325100438 0.351 0.382 — — — — —

Fig. 3. The smeared interface between oil and water.

would fall in the given capillary number range and to also ensure that the experimental values at those capillary numbers are available in the literature. The studied flow geometry is shown in Fig. 1.

The simulation results are summarized in Table 1. The experimental results, originally taken from four different papers, have been tabulated by Ref. [5] into giant table. Table 1 shows the lower and upper bound film thickness obtained from the simulations for each capillary number and compared them to the published experimental results and to Equation (5) for small capillary numbers. Different papers give film thicknesses for different values of capillary numbers and occasionally more than one paper has reported a film thickness for the same capillary number. Entries with dashes indicate unavailable data for that capillary number. Table 1 results are depicted in a graphical form in Fig. 4. According to Fig. 4, simulation and experimental results are quite close to each other.

Experimental values from two different researches can differ somewhat, perhaps due to measurement errors. This difference can be depicted by the purple and orange dots representing reference [13] and [5], respectively in Fig. 4. Our simulation

results appear to match more closely to reference [13] data than reference [5].

The authors of Ref. [5] have proposed from their theory that as capillary numbers become larger, the film thickness will level off at a constant value of 5 = 0.2. Their experiments have found a leveled off value of approximately 0.27. The author of Ref. [13] found the leveled off value to be approximately 0.32. According to the simulation herein, this value seems to be about 0.36 and the leveling off phenomena can be seen in Fig. 4.

3.4. Reynolds number effect

The small discrepancy between our simulation and published experimental values can be explained for by considering the Reynolds number effect. We can see from equation (2) that a same capillary number can be derived with a different combination of input parameter on the right hand side of the equation. Those same input parameters are included in the Reynolds number calculation as shown in equation (3). Hence, same capillary number can have different Reynolds number. It has been shown by Refs. [9] and [11] that Reynolds number can have an effect on the film thickness and since we have only matched

Table 2

Effect of Reynolds number on the film thickness at constant capillary number and viscosity ratio.

Ca Re Pwater (kg/m3) Poil (kg/m3) mwater (kg/(m S)) moil (kg/(m S)) Yoil-water (N/m) M Vinlet (m/s) AFT (5)

Case 1

Simulation 1 0.00096 49.56 998.2 860 0.001003 0.001003 0.1215 1.0 0.1156015814 0.0160

Simulation 2 0.00096 1.74 998.2 860 0.001003 0.001003 0.003375 1.0 0.0040586512 0.0164

Case 2

Simulation 1 0.01 14.43 998.2 860 0.001003 0.001003 0.003375 1.0 0.033658814 0.0577

Simulation 2 0.01 28.85 998.2 860 0.001003 0.001003 0.00675 1.0 0.0672943023 0.0582

Simulation 3 0.01 64.93 998.2 860 0.001003 0.001003 0.01519 1.0 0.151453 0.0588

Table 3

Effect of viscosity ratio on the film thickness at a fixed capillary number and low Reynolds number.

Ca M Pwater (kg/m3) Poil (kg/m3) mwater (kg/(m S)) moil (kg/(m S)) Yoil-water (N/m) Vinlet (m/s) Re AFT (5)

Case 1

Simulation 1 0.2 39.9 998.2 860 0.001003 0.04 0.1215 0.1348837209 1.45 0.2524

Simulation 2 0.2 73.8 998.2 860 0.001003 0.074 0.003375 0.1824186047 1.06 0.2519

Simulation 3 0.2 147.6 998.2 860 0.001003 0.148 0.10125 0.1376744186 0.40 0.2521

Case 2

Simulation 1 2 39.9 998.2 860 0.001003 0.04 0.003375 0.168372093 1.81 0.3663

Simulation 2 2 73.8 998.2 860 0.001003 0.074 0.00675 0.1824186047 1.06 0.3575

Simulation 3 2 99.7 998.2 860 0.001003 0.1 0.01519 0.1348837209 0.58 0.3561

the capillary numbers and not the Reynolds number, it could very well be the reason for the small discrepancy. This is perhaps the reason why there is a larger error between our simulation results and that of Ref. [5]. The authors of Ref. [5] have not reported the velocity or flow rate used in the experiment. Hence, it is not possible to calculate their Reynolds number and compare it with ours. Fig. 2 in Ref. [11] and Fig. 10 in Ref. [9] shows that for capillary number less than or equal to 0.01, the Reynolds number has very little effect in the range of 0—70 and the film thickness increases with increasing Reynolds number but the increase is very small and it almost looks as if the film thickness is independent of the Reynolds number.

To test if this is indeed the case and to see the Reynolds number effect, two different simulations were run. Water density and viscosity are kept constant, as shown in Table 2 and oil properties and interfacial tension are changed to generate different capillary and Reynolds numbers. In Table 2, the results of two different cases involving the same capillary number, at different Reynolds numbers, are shown. Case 1 is for capillary number 0.00096 and Case 2 is for capillary number 0.01. The viscosity ratio for all simulations is kept constant at 1 to mitigate its effect. According to the table, it can be observed that for any given capillary number, the film thickness becomes thicker as the Reynolds number is increased but the increase is very small. This is exactly the prediction of the authors of Refs. [9,11 ]. Although Refs. [9,11 ] have obtained their results as for a gas bubble displacing a liquid, according to the above table, a similar trend seems plausible for the flow of liquid displacing another liquid. This trend would be confirmed if similar simulations are ran at

i 0.350

: 0.300

; 0.250

: 0.200 J

I 0.150 j 0.100 ' 0.050 0.000

1.0 1.5

Capillary Number

♦ Simulation lower bound ■ Simulation upper bound Aul and Olbricht • Carvahlo et al.

• Goldsmith and Mason • Beresnev et al. Bretherton's Equation

Fig. 4. Ca vs. normalized film thickness. Comparison of the film thickness values from the simulation with that of published experimental values for the capillary number range of 0.00096-2.

higher capillary numbers but due to the time constraints as well as high inlet velocities at higher capillary numbers leading to unstable VOF method simulations, this was not possible.

3.5. Viscosity ratio effect

According to Ref. [5], there is no noticeable change in the film thickness with changing Viscosity ratio. However Ref. [25], predicts thinning of the film thickness with the increase in the Viscosity ratio. The conflicting piece of information in the two papers is the disagreement between the effects of viscosity ratio on wetting film thickness. Hence, we would like to predict how Viscosity ratio affects the film thickness from our simu-lations.The authors of Ref. [25] have said that the thinning of film thickness with increase of viscosity ratio become more apparent at Ca > 1 based on the simulations. At lower capillary numbers, the effect is still there but it hard to detect because of the experimental uncertainty and could be the reason why it is not seen by some authors. They further state that at high capillary number experiments cannot be carried out because of the limitation of the experimental procedures.

We have run three different simulations for capillary number 0.2 and 2 at different viscosity ratio and kept the Reynolds number less than two to minimize the inertial effects. The Reynolds number cannot be kept constant because it also involves the changing parameter of viscosity. Table 3 shows the parameters used in the three simulations. The highest viscosity ratio we can use is about 148. Simulations employing viscosity ratios higher than 148 produce incorrect results and it is one of the shortcomings of the VOF method. From the table, we can see that the film thicknesses for all three cases of capillary number 0.2 are extremely close to each other and we cannot see any definitive trend in the film thickness with increasing viscosity ratio. For the capillary number of 2 we notice a similar phenomenon. Our results seem to agree with the conclusion of [5] that the film thickness is independent of the Viscosity ratio. Simulations at capillary numbers in between 0.2 and 2 and higher viscosity ratio must be conducted to see if results of ref. [5] are valid for a large range of capillary numbers and viscosity ratio. However, due to time constraints and shortcoming of VOF method at higher viscosity ratio than 148, this was not possible.

4. Conclusions

An immiscible liquid—liquid multiphase flow simulations using volume of fluid method have been conducted. For the range of capillary numbers studied, a good match between the simulation and published literature values are found. For the range of viscosity ratios studied, the film thickness seems to be independent of the viscosity ratio. The effect of Reynolds number on the film thickness as described by Refs. [9,11] also seems plausible for liquid—liquid flow from our simulations. A good

understanding of a pore—scale water flooding mechanism is obtained.

Volume of fluid method can have some potential factors that are worth considering for the future simulations. One of the major concern for the future work is that the properties of two fluids used in the multiphase flow such as density and viscosity cannot be very different from each other [16]. These can limit the types of simulations that can be performed. For example Ref. [5], has performed core-annular flow experiments with viscosity ratios of 3650. To perform simulations with high viscosity ratio with VOF method would be difficult due to the large viscosity differences between the two liquids. Due to time constraints, only limited viscosity ratios have been tried and simulations with higher viscosity ratios should be tried in the future work. The capillary numbers in the simulations are adjusted based on changing velocity, viscosity, and the interfacial tension. High velocity lead to high courant number and this in turn lead to stability issues that require small time steps. As mentioned before, due to the nature of the VOF method, the ranges of the liquid properties that can be employed in the simulation need to be explored so that the maximum capillary number and viscosity ratio that can be simulated is known.

The core—annular flow simulations can generate reliable results and if used correctly, they can save a lot of time and effort needed to perform equivalent experiments. The mesh used in the numerical development in this paper will be used further to simulate single oil slug flow surrounded by water. These simulations would be representative of a situation of left over oil after water flooding.

Acknowledgments

The authors are grateful to International Performance Assessment Centre for Geologic Storage of Carbon Dioxide (IPAC—CO2) for allowing us to use that high performance cluster for the computations. We would also like to thank Patrick Mann of IPAC—CO2 for his help with the setup of Fluent on the cluster and for his continual help with the remote operation of Fluent.

Nomenclature

Ca capillary number

Re Reynolds number

Co courant number

Ax grid size, m

At time step, s

Vfluid fluid velocity, m/s

M Viscosity ratio

moil viscosity of oil, kg/m-s

mwater viscosity of water, kg/m-s

mannular viscosity of the fluid occupying annuals of the tube, kg/ m- s

mcore viscosity of the fluid occupying core of the tube, kg/ m- s

poil density of oil, kg/m3

D diameter of the tube, m

5 normalized film thickness

Vinlet inlet velocity, m/s

goil-water oil—water interfacial tension, N/m

UB upper bound

LB lower bound

AFT average film thickness = (LB + UB)/2

References

[1] F. Fairbrother, A.E. Stubbs, The bubble-tube methods of measurement, J. Colloid Sci. 1 (0) (1935) 527—529.

[2] F.P. Bretherton, The motion of long bubbles in tubes, J. Fluid Mech. 10 (2) (1961) 166—188.

[3] L.W. Schwartz, H.M. Princen, A.D. Kiss, On the motion of bubbles in capillary tubes, J. Fluid Mech. 172 (1986) 259—275.

[4] J. Ratulowski, H.-C. Chang, Marangoni effects of trace impurities on the motion of long gas bubbles in capillaries, J. Fluid Mech. 210 (1990) 303—328.

[5] I. Beresnev, W. Gaul, R. Dennis Vigil, Thickness of residual wetting film in liquid-liquid displacement, Phys. Rev. E. 84 (2) (2011).

[6] R.N. Marchessault, S.G. Mason, Flow of entrapped bubbles through a capillary, Ind. Eng. Chem. Res. 52 (1) (1960) 79—84.

[7] G.I. Taylor, Deposition of a viscous fluid on the wall of a tube, J. Fluid Mech. 10 (2) (1961) 161—165.

[8] J. Chen, Measuring the film thickness surrounding a bubble inside a capillary, J. Colloid Interface Sci. 109 (2) (1986) 341—349.

[9] M.D. Giavedoni, F.A. Saita, The axisymmetric and plane cases of a gas phase steadily displacing a newtonian liquida simultaneous solution of the governing equations, Phys. Fluids 9 (8) (1997) 2420—2429.

[10] S. Irandoust, B. Anderson, Liquid film in taylor flow through a capillary, Ind. Eng. Chem. Res. 28 (11) (1989) 1684—1688.

[11] M. Heil, Finite reynolds number effects in the bretherton problem, Phys. Fluids 13 (9) (2001) 1—6.

[12] R.W. Aul, W.L. Olbricht, Stability of a thin annular film in pressure-driven, low-reynolds-number flow thorough a capillary, J. Fluid Mech. 215 (1990) 585—599.

[13] E.J. Soares, M.S. Carvahlo, P.R. Souza Mendes, Immiscible liquid-liquid displacement in capillary tubes, J. Fluids Eng. 127 (2005) 24—27.

[14] H.L. Goldsmith, S.G. Mason, The flow of suspensions through tubes. ii single large tubes, J. Colloid Sci. 18 (3) (1963) 237—261.

[15] S.R. Hodges, O.E. Jensen, J.M. Rallison, The motion of a viscous drop through a cylindrical tube, J. Fluid Mech. 501 (2004) 279—301.

[16] G. Wang, Wave Scattering and Propagation in Porous Media with Excitations from Multiple Wave Sources, University of Regina, Regina, Canada, 2008 [PhD thesis].

[17] G. Cheng, Experimental and Numerical Study of Residual Oil Mobilization Under Static and Dynamic Conditions, University of Regina, Regina, Canada, 2009 [Master's thesis].

[18] M. Renardy, Y. Renardy, J. Li, Numerical simulation of moving contact line problems using a volume-of-fluid method, J. Comput. Phys. 171 (1) (2001) 243—263.

[19] C.W. Hirt, B.D. Nicholas, Volume of fluid (vof) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1) (1981) 201—225.

[20] Fluent 13.0 User's Guide, 2010. Technical report.

[21] R. Gupta, D.F. Fletcher, B.S. Haynes, On the cfd modeling of taylor flow in microchannels, Chem. Eng. Sci. 64 (12) (2009) 2941 —2950.

[22] D.L. Youngs, Time-dependent multi-material flow with large fluid distortion, in: K.W. Morton, M.J. Baines (Eds.), Numerical Methods for Fluid Dynamics, 1982, pp. 273—285.

[23] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (2) (1992) 335—354.

[24] A. Bayliss, C.I. Goldstein, E. Turkel, A method of determining the thickness of liquid—liquid interfaces, Colloids Surfaces 113 (1—2) (1996) 51—59.

[25] E.J. Soares, R.L. Thompson, Flow regimes for the immiscible liquid liquid displacement in capillary tubes with complete wetting of the displaced liquid, J. Fluid Mech. 641 (2009) 63—84.

[26] H. Joshi, A Numerical Approach for Oil Slug Mobilization in a Capillary Tube under Absence and Presence of External Excitations, University of Regina, Regina, Canada, 2013 [Master's thesis].