Computational Fluid Dynamics to develop novel correlations for residual saturation of the displaced fluid in a capillary tube

Journal Pre-proof Computational Fluid Dynamics to develop novel correlations for residual saturation of the displaced fluid in a capillary tube

Meisam Adibifard, Ali Nabizadeh, Mohammad Sharifi PII:

S0167-7322(19)33622-0

DOI:

https://doi.org/10.1016/j.molliq.2019.112122

Reference:

MOLLIQ 112122

To appear in:

Journal of Molecular Liquids

Received date:

29 June 2019

Revised date:

6 October 2019

Accepted date:

11 November 2019

Please cite this article as: M. Adibifard, A. Nabizadeh and M. Sharifi, Computational Fluid Dynamics to develop novel correlations for residual saturation of the displaced fluid in a capillary tube, Journal of Molecular Liquids(2018), https://doi.org/10.1016/ j.molliq.2019.112122

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2018 Published by Elsevier.

Journal Pre-proof Computational Fluid Dynamics to Develop Novel Correlations for Residual Saturation of the Displaced Fluid in a Capillary Tube

1 2

Meisam Adibifarda*, Ali Nabizadehb, Mohammad Sharific

3

a

4

Swalm School of Chemical Engineering, Mississippi State University, MS, United States

b

Institute of Petroleum Engineering, School of Energy, Geoscience, Infrastructure and Society, Heriot-Watt University, Edinburgh, UK

5 6

c

7

Petroleum Engineering Department, Amirkabir University of Technology, Tehran, Iran

8

oo

f

*E-mail: [email protected]

pr

Abstract:

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

Keywords: Computational Fluid Dynamics; Finite Volume; Volume of Fluid; Immiscible Displacement; Capillary Flow; Dynamic Contact Angle

32 33

Jo u

rn

al

Pr

e-

The main goal of this study is to numerically investigate the immiscible fluid-fluid displacement at the micro-level under low inertial and buoyancy forces. To this end, a micro capillary tube was constructed, and the fluid-fluid interface and fluid-fluid-wall triple line were tracked using VOF (Volume of Fluid)- CSF (Continuum Surface Force) strategy. Meshindependency analysis was conducted by utilizing the Adaptive Mesh-Refinement (AMR) technique. Simulation results agreed well with the experimental data obtained from literature with the AARD (Average Absolute Relative Deviation) of around 10%. The validated model was then used to generate the CDC (Capillary Desaturation Curve) curves for different wall wettabilities, considering the dynamic contact angle in the simulations. The generated CDC curves at the breakthrough time showed that there are three different zones, namely low-Ca, transitional, and high-Ca region. The only difference between saturation profiles of different wettabilities was in the transition zone, based on the CDC analysis. The saturation and velocity profiles at the central plane of the tube were constructed and compared. Phase separation can happen depending on Ca number and the magnitude of the dynamic contact angle. Eventually, outcomes of the numerical simulator were used to develop new correlations for the residual saturation as functions of dynamic contact angle and Ca number. Nonlinear regression analysis with Levenberg-Marquart algorithm was used to find coefficients of the correlations. The developed correlations showed to be highly robust for the transition and high-Ca regions with the correlation coefficients (R2) up to 0.96 and the maximum AARD of 12%. Results of this study are applicable to both capillary-fingering and viscous-fingering flow-regimes based on the Lenormand‟s phase-diagrams.

34

Journal Pre-proof 35

Various industrial processes rely upon stratified fluid-fluid displacements [1].

36

„Coextrusion‟ in the polymer industry [2], multi-phase flow within the soil [3], cementation

37

of oil reservoirs [4, 5], and Enhanced Oil Recovery (EOR) in the petroleum industry [4, 6]

38

are few instances in which immiscible fluid-fluid displacement is of great importance. Precise

39

characterization of the immiscible fluid displacement is crucial in the aforementioned

40

processes because it directly influences the efficiency of the process.

41

oo

f

1. Introduction

42

chief factor behind the remained oil saturation after water-oil displacements [7]. To decrease

43

pr

It is well known that capillary trapping of the non-wetting phase by wetting phase is the

e-

the saturation of the trapped oil by increasing the Ca (Capillary) number, several EOR

Pr

schemes such as surfactant injection and ASP (Alkali-Surfactant-Polymer) flooding have

44 45 46

immiscible micro-flow displacements; therefore, unraveling the underlying relationship

47

al

been proposed in the literature [8, 9]. In other words, Ca number plays a substantial role in

rn

between the volume of the trapped fluid and Ca number will assist engineers in different

Jo u

fields to rigorously design their own immiscible fluid-fluid experiments.

48 49

Numerous experimental and numerical models have been utilized in the literature to

50

qualify two-phase flow displacement in capillary tubes and other geometrical conducts. For

51

example, Goldsmith and Mason (1963) investigated liquid-liquid immiscible displacement

52

under laminar flow conditions within circular tubes. They considered effects of Ca number

53

and viscosity ratio on the magnitude of the mass of liquid attached to the wall [10]. A micro-

54

scale experimental study using X-ray micro-tomography was conducted by Iglauer et al.

55

(2012) to monitor residual oil saturation in sandstone reservoirs [11]. An experimental

56

research was proposed by Lenormand et al. (1983) to analyze immiscible two-phase flow

57

within capillaries. In their experiments, they found films of wetting fluids in the corners of

58

Journal Pre-proof 59

methods with experimental approaches to investigate immiscible drainage flow in

60

micromodels. They proposed three domains for the drainage displacement including capillary

61

fingering, viscous fingering, and stable displacement [13]. Likewise, Zendehboudi et al.

62

(2012) proposed a coupled experimental-computational approach to simulate gravity drainage

63

processes in naturally fractured reservoirs. The impacts of crucial parameters, including

64

fractures‟ aperture, length, and orientation, were investigated in their study [14].

65

f

the capillaries [12]. In another study, Lenormand et al. (1988) coupled computational

66

and sandstones to probe the effects of pore structure and heterogeneity over saturation of the

67

pr

oo

Chatzis et al., (1983) carried out several tests using glass bead packs, 2D-micromodels,

68

scheme to simulate the flow of inviscid bubbles and viscous drops in capillary tubes [16]. Liu

69

and Gao (2007) investigated the two-phase flow of water-gas in vertical and horizontal

70

capillary tubes. Their results showed that surface tension is of great importance for wavy

71

stratified flow regimes [17]. Feng (2009) used Galerkin Finite-Element technique to examine

72

the steady axisymmetric behavior of a long gas bubble moving within a liquid in a capillary

73

tube [18]. Freitas et al. (2011) used the same numerical approach utilized by Feng (2009) to

74

Jo u

rn

al

Pr

e-

trapped phase [15]. Westborg and Hassager (1989) employed Finite-Element numerical

analyze the flow of two immiscible liquids in a plane channel. They ended up with the

75

conclusion that under low Ca number conditions, the obtained analytical results for gas-

76

displacement process in a capillary tube can be used for the plane channel as well [4].

77

Last year, Adibifard (2018) introduced a new analytical solution that can be used to

78

reconstruct the saturation-time data solely by matching the experimental inlet pressure data.

79

The developed techniques proved to be robust with maximum relative error of around 11%

80

[19]. Recently, Nabizadeh et al. (2019) studied the effects of dynamic contact angle on

81

multiphase displacement flow in angular pores and showed that implementing dynamic

82

Journal Pre-proof contact angle models can significantly decrease the discrepancy between computational

83

analysis and experimental observations [20].

84 85

Stokes equation can be used accurately to characterize the fluid flow field rather than the

86

renowned Darcy law [21] which is usually used to simulate flow in porous media. Several

87

research studies have utilized pore-scale modeling to simulate the multi-phase flow in micro-

88

scaled porous media. Kovscek et al. (1993) proposed a pore-scale investigation to simulate

89

f

It is well established that at the small length scale, i.e. scale of a single pore, Navier-

90

modeling study to investigate the relationship between empirical wettability measurement

91

pr

oo

multiphase flow through porous media [22]. Dixit et al. (2000) developed a pore-scale

92

of wettability on multiphase flow [24].

93

Pr

e-

methods [23]. Blunt (1998) proposed another pore-scale study to demonstrate the importance

94

equation can be used to precisely describe flow in different scales and different geometries.

95

al

CFD (Computational Fluid Dynamics) is a field of engineering in which Navier-Stokes

rn

CFD has gained much attentions in investigating the different displacement processes under

96 97

simulation of gas-liquid-solid three phase flow (Takahashi et al., 2016), simulation of soil-

98

erosion through coupled CFD-DEM (Discrete Element Method) [26], examining performance

99

of isothermal and non-isothermal EOR schemes [27], and analyzing the effects of initial

100

wetting film on immiscible flow in pore-doublet [28].

101

Jo u

different fluid flow circumstances. Examples are studying elasto-visco-plastic fluids [25],

Given the importance of CFD modeling in micro-flow simulations, in this study, we use

102

CFD approach to simulate the immiscible multi-phase flow in a micro capillary tube

103

subjected to certain boundary conditions. The main goal of this study is to investigate the

104

effects of fluid-fluid-wall contact angle and Ca number on the residual saturation of the

105

displaced fluid, after model validations using experimental data. To capture the fluid

106

Journal Pre-proof dynamics in the tube, we consider dynamic contact angle in our calculations. Plotting CDC

107

(Capillary Desaturation Curve) curves for different wettabilities, we propose three different

108

correlations to estimate the residual saturation of the non-wetting fluid. The proposed

109

correlations are very robust for the examined Ca region and can be used to regenerate CDC

110

curves using Ca number and tube‟s wettability, where there is lack of enough experimental

111

data.

112

oo

f

113

2. The Developed Model

114 115

single capillary tube. A commercial CFD software was employed to construct the geometry,

116

e-

pr

A FV (Finite Volume) numerical scheme is used in this study to simulate fluid flow in a

Pr

discretize and apply the boundary conditions, and run the FV simulations. Although there are

117 118

Element), and XFE (Extended Finite Element) [29-32], the chief advantages of the FV

119

al

many numerical schemes proposed in the literature such as FD (Finite Difference), FE (Finite

rn

technique over other numerical schemes are [33]: Its applicability for unstructured gridding



Using integral formulation of conservation law

Jo u



120 121 122

Simulation of an immiscible displacement requires solving for a free-boundary problem.

123

Different techniques have been proposed in the literature for this purpose, and each of them

124

can be used for a specific case. In general, the methods employed to locate the fluid-fluid

125

interface include VOF (Volume of Fluid), ASM (Algebraic Slip Mixture), Eulerian, and

126

Lagrangian [34].

127

VOF technique lies within the tracking interface methods which are used to describe the

128

behavior of multi-phase flow by considering a single set of Navier-Stokes equation for the

129

Journal Pre-proof 130

technique, fluid properties of the total system are calculated by adopting a simple averaging

131

method on the volume fraction and the properties of the phases. On the other hand, in the

132

Eulerian technique, Navier-Stokes equation should be rewritten for each flowing phase. As

133

described by [34], the Eulerian method is primarily favorable to model the liquid-liquid

134

dispersion phenomena in contactors. Also, the Lagrangian technique is mainly appropriate for

135

the particles‟ transfer through a continuous phase as this method takes into account the very

136

small volume fractions for the dispersed phase [34].

137

oo

f

total system and track the interface by assigning a volume fraction to each cell. In this

pr

In this study, we utilize VOF technique to track the oil-water free interface and CSF

138 139

utilize an empirical equation proposed in the literature to take into account the dynamic

140

contact angle in our simulations. The dynamic contact angle is calculated using Ca number

141

and is used within the VOF-CSF technique in our numerical simulations. In the subsequent

142

subsections, we elaborate on the mathematical formulations of our VOF-CSF scheme and

143

dynamic contact angle calculations, as well as the transient modeling of the simulation.

144

Jo u

rn

al

Pr

e-

(Continuum Surface Force) method to consider the fluid-fluid-wall contact angle. We also

145

2.1. VOF-CSF Strategy

146

VOF tracking strategy has been broadly used in the literature to simulate multi-phase

147

flow in different industrial processes such as fuel cells [35, 36], centrifugal contactors [37,

148

38], multi-phase flow in pipe [39, 40]; ligament-growth [41]; cavitating flow [42, 43], and

149

heat and mass transfer modeling in a deformable interface [44]. VOF is an effective and

150

adaptable method to consider the effects of free-surface in transient flow analysis, even in

151

dealing with complicated free-boundary problems [45, 46]. Another important advantage of

152

Journal Pre-proof VOF is that it assures the precise mass conversation and can be used for any mesh geometry,

153

in addition [47]. The mathematical formula for this technique is provided in the following.

154

Following is the main momentum equation used in the VOF technique to determine the

156

velocity field and pressure for the total multi-phase system [48]: (

) ̿

( )

⃗

155

(1)

157

158

f

(2)

159

stress tensor, t is time, is the fluid‟s density, and u is the fluid‟s velocity. The first equation is

160

pr

oo

Where p is the local pressure, F(σ) is the capillary force imposed to the interface, is viscous

e-

the momentum or Navier-Stokes equation while the second one is the continuity equation for

161 162

for both continuity and momentum equations.

163

Pr

an incompressible fluid. The scaled residual function was used as the convergence criterion

al

Then, the color function, or volume fraction, of the system is calculated using the

165

(3)

Jo u

rn

following transport equation [48]:

⃗

164

Where stands for the color function representing the local volume fraction of the first phase in the VOF technique.

166

167 168

In addition to the VOF strategy, Continuum Surface Force (CSF) is used to incorporate

169

the interfacial forces between the two-phases as a source term in the general momentum

170

equation, as seen before in equation (1). This method was developed by [49] based on the

171

continuum method which eliminates the interfacial reconstruction [49]. Since constant

172

Interfacial Tension (IFT) is assumed to model the interface‟s curvature, CSF model can be a

173

good option for our simulations [50]. Surface‟s wettability was considered through imposing

174

Journal Pre-proof a contact angle between the fluid-fluid interface and the solid wall. The fluid-solid contact

175

angle is treated as a dynamic boundary condition and influences the curvature of the interface

176

in the cells juxtaposed to the solid wall [50].

177

Capillary force F(σ) in the Navier-Stokes equation is calculated using the following equation:

178

, ( )

(

)

179

(4)

180

and n is the normal vector of the interface, defined as following:

181

oo

f

Where  refers to the Dirac delta function, is the interfacial tension, is interface‟s curvature,

(5)

pr

,,

e-

Where is the color functuion calculated from the PDE (Partial Differential Equation) in

182 183 184

calculated using the following equation:

185

rn

al

Pr

equation (3). The interface‟s curvature  should be calculated too. This parameter is

| | ,

186

(6)

187

flow simulations. The fluid-fluid-wall triple line is taken into consideration by adjusting the

188

normal vectyor nearby the wall. Therefore, following equation is used to update the surface

189

normal vectors for cells next to the wall:

190

Jo u

Wall‟s wettablity or contact angle is another parameter that significantly impacts the fluid







n  n w cos   t w sin  ,

(7)

is the unit vector normal to the wall and is the unit vector tangential to the wall.

191 192

A geometric reconstruction method is used within the VOF technique to find the fluid

193

fluxes at the cells adjacent to the fluid-fluid interface. A linear piecewise approximation is

194

Journal Pre-proof used to represent the two-phase interface, and then the linear slope at each cell is used to find

195

the advection of fluids through the cell‟s faces [50].

196 197 198

Using VOF-CSF scheme, we can only enter one value for the contact angle of fluid-fluid-

199

wall. As it is understood, contact angle of the fluid-wall under dynamic conditions is different

200

from the contact angle measured under static conditions in the laboratory. Hence, to consider

201

oo

f

2.2. Dynamic Contact Angle

202

based on static contact angles and Ca number by using the equation adopted from [51]. Then,

203

the calculated dynamic contact angle was used in the CFD simulator for our simulations. The

204

e-

pr

effects of the fluid‟s dynamic on the contact angle, we calculated the dynamic contact angle

(

)

( )

)

(

),

,

(8)

al

(

Pr

equation used in this study to calculate dynamic contact angle is as following[51]:

rn

Where and are the static advancing and dynamic advancing contact angles, respectively.

Jo u

Also, Ca is the capillary number governing the flow regime.

205

206

207 208

Using Ca number and static contact angles used in this study (i.e. 0, 90, and 180), the

209

corresponding contact angles for the dynamic condition were calculated. In this method, we

210

assumed an average dynamic contact angle, and it might be considered as an approximated

211

dynamic contact angle analysis. We assume that the length of the capillary tube is short

212

enough to make such approximation. As it is obvious from equation (8), the dynamic and

213

static contact angles will remain the same for the capillary tube with = 180, i.e. oil-wetting

214

medium.

215 216

2.3. Transient Modeling

217

Journal Pre-proof An Explicit scheme was used to update the solution in time; the magnitude of the time-

218

step was adjusted automatically during the simulations by continuously monitoring the

219

Courant number within the solution domain. The maximum Courant number of two was used

220

as a cutoff to restrict the time-step and also to avoid large magnitudes of inter-cell fluid

221

velocities. The Courant number is part of the CFL (Courant-Friedrichs-Lewy) condition

222

introduced by [52] and is defined as flowing for any arbitrary n-dimensional problem:

223

⃗

,,

(9)

224

oo

f

⃗⃗⃗⃗⃗

225

space. Applying the CFL conditions, the time increment Δt is calculated using the following

226

equation [50]:

227

e-

pr

Where is the fluid velocity, and stands for the spatial increments in the three dimensional

)

Pr

(

,

(10)

al

Where is the spatial increment and is the maximum of local eigenvalues of a preconditioned

rn

system used to modify the time-derivative of the governing equations [50]. Finally,

Jo u

Incompressible flow conditions were assumed for both displacing and displaced fluids because of the negligible density changes in micro-scale.

228

229 230 231 232

The graphic presented in Fig. 1 at best depicts the major steps used in this study to

233

develop our novel correlations based on the numerical results obtained from CFD

234

simulations. These steps are followed in the subsequent sections within the text.

235 236 237

Figure 1. A schematic of all the main steps utilized in this study to develop new

238

correlations for residual saturation of the displaced fluid

239

Journal Pre-proof 240 241 242 243

In the first step, mesh-independence analysis was conducted to show that the system‟s

244

solution is not affected by the number of nodes and/or elements. Subsequently, to show that

245

f

3. Results and Discussion

246

regenerated using the numerical model. Two-phase flow simulations were repeated then

247

pr

oo

results are physically reliable, experimental data belonging to the water-oil displacement was

248

curves were generated for different fluid-wall wettabilities at the time of breakthrough.

249

Finally, three correlations were developed to estimate the residual oil saturation from Ca

250

Pr

e-

under different Ca number and fluid-wall contact angles. Based on the obtained results, CDC

251

and oil-wetting mediums.

252

Jo u

rn

al

number and dynamic contact angle, for three different static contact angles, i.e. water, neutral

253

3.1. The Geometry and Boundary Conditions

254

In this study, we used a capillary tube with the length of 200 μm and radius of 5 μm. The

255

discretized geometry is illustrated in Fig. 2. To reduce the number of cells, the sweeping

256

technique was used to sweep the triangular meshes generated in the 2D inlet surface.

257

Accordingly, unstructured prism cells were used to discretize the fluid domain.

258

Figure 2. The meshed capillary tube used for the numerical FV simulations

259 260

Journal Pre-proof Constant inlet velocity was applied as a boundary condition to the inlet face whist

261

pressure outlet of one atmosphere, i.e. 14.7 psi, was imposed to the outlet face of the tube. A

262

no-flow boundary condition was imposed over the peripheral surface area of the tube, which

263

implies no slipping fluid at the fluid-solid contact area.

264 265

defined by [56], which the most applicable definitions of Ca number [60], was used in this

266

study. Therefore, Ca number and Nμ are defined as following:

267

(11)

268

(12)

pr

269

e-

oo

f

Although there are several definitions for Ca number in the literature [53-59], the one

270

Pr

In the above equations, indices 1 and 2 refer, respectively, to the displacing and displaced

271

IFT (Interfacial Tension) acting between two fluids. Capillary number represents the relative

272

al

fluids. is the viscosity term; is the average front velocity for the displacing fluid, and is the

rn

importance of the viscous or capillary forces in the conduit [61, 62]; Lower Ca numbers

273 274

represent flows with viscous forces as the driving force.

275

Jo u

indicates that capillary is the main driving force for the fluid flow while larger Ca numbers

The properties of the injection fluid are adopted from [63] and are listed in Table 1. The

276

oil phase has viscosity of 0.05233 Pa.s (52.33 cp) and density of 920 kg/m3 (Soares et al.,

277

2015).

278 Table 1. Properties of the injection fluids used in this study adopted from [63]

Fluid

Nµ

IFT, mN/m

F1C

8.99

5.30

279

280

3.2. Mesh Independence Analysis

281

Journal Pre-proof 282

the solution from the grid size. In fact, matching numerical results with experimental data

283

does not necessarily mean that the model‟s output is independent from the discretized

284

geometry [64]. Given that, performing mesh sensitivity analyses to obtain a grid-independent

285

model is a critical step for any CFD simulation. Although various approaches have been

286

suggested to perform grid independence analysis [64], the AMR (Adaptive Mesh

287

Refinement) technique, which is a dynamic method introduced by Berger and Colella (1989)

288

[65], is employed in the current study.

289

oo

f

Mesh independence analysis is ineluctable in order to demonstrate the independency of

pr

To carry out the mesh independency analysis, simulation was conducted at Ca≈0.1.

290 291

displacement process. After solving for the velocity and pressure fields using contact angle of

292

1 0 for the water phase, i.e. oil wetting medium, the second spatial derivative of the pressure

293

at each single cell was used to adjust the number of cells in the model. Using this procedure,

294

dimensions of the cells with higher pressure curvatures were shrunken to precisely capture

295

the pressure profile within the domain.

296

rn

al

Pr

e-

Initially, the coarse meshing shown in Fig. 2 was used to simulate the fluid-fluid

297

initial value of 14553 to 82278, which resulted in huge increase in the computational time.

298

Although the number of cells was increased substantially, the average inlet pressure and the

299

oil saturation data plotted in Fig. 3 showed trivial differences between the fine and coarse

300

mesh scenarios, particularly for the times before the breakthrough. As a result, the coarse

301

meshing was accurate enough for fluid flow simulations and therefore was used for the

302

subsequent simulations. The important properties of the meshes such as orthogonal quality

303

and aspect ratio are provided in Table 2.

304

Jo u

Conducting the mesh refinement process, the number of cells were increased from its

305

Journal Pre-proof Figure 3. Average pressure at the inlet face, and the average oil saturation within the capillary tube versus time for coarse and adapted meshing geometries

306 307 308 309 310 311

oo

f

312

Mesh Statistic

Maximum aspect ratio

Pr

Minimum cell volume

e-

Minimum orthogonal quality

pr

Table 2. Mesh statistics for the coarse meshing

0.661 5.932 5.809e-19 m3 1.700e-18 m3

rn

3.2. Model Verification

Value

al

Maximum cell volume

313

314 315 316

the open literature, were regenerated by using the developed model. To this end, two-phase

317

flow data belonging to a straight capillary tube was collected from [1], and the residual oil

318

saturation data at different Ca number were used for the verification purpose. It is worth to

319

mention that the actual geometry of the experimental model was transformed to a micro-level

320

by preserving the aspect ratio of the original system (see Fig. 2). Since Ca number solely

321

depends upon the physical and chemical properties of the involved fluids, its magnitude is not

322

affected by the corresponding scale transformation. Additionally, since low Re number

323

condition is maintained for the reduced size model, it is expected to observe the same

324

displacement patterns for the scaled model.

325

Jo u

To verify the accuracy of the numerical model, a set of experimental data, collected from

Journal Pre-proof Simulations were conducted for different Ca number at viscosity ratio of 8.99. Since [1]

326

measured the residual oil saturations at the 7% remained to the end of the tube, we measured

327

the residual volume fraction at distance 1 6 μm. Since original experimental data was

328

conducted in an oil-wetting tube, the dynamic and static contact angles would be the same in

329

this case as discussed before in the modeling section.

330 331

obtained between the experimental and simulation results with Average Absolute Relative

332

f

Simulation results are compared with the experimental data in Fig. 4. A good match was

333

experimental data, the numerical model can accurately predict the trend in the saturation data

334

pr

oo

Deviation (AARD) of 10.27%. Despite the existing deviations between the numerical and the

335

reaches to a stabilized value for the large values of Ca number.

336

al

Pr

e-

against time; numerical residual oil saturation increases by increasing the Ca number and

Jo u

rn

Figure 4. A comparison between the experimental and FV numerical results based on the saturation-Ca curves

337 338 339 340 341

Detailed statistical analysis of the numerical results is provided in Table 3. According to

342

this analysis, the correlation coefficient (R2) is as high as 0.93 and the error standard

343

deviation is as low as than 0.0356, which shows a good match between the experimental and

344

numerical data.

345 346

Ultimately, graphical match between the numerical and experimental data, accuracy of

347

the numerical model in predicting the trend of the saturation-Ca curve, high values of the R-

348

squared, and low values of the STD demonstrate the robustness and reliability of the

349

Journal Pre-proof developed model in prediction of the multi-phase saturation map. The verified model is used

350

in the next section to generate the CDC curves for different fluid-solid equilibrium lines.

351 352

Table 3. Detailed Statistical Analysis of the numerical outcomes

Injection Fluid 8.99

0.9356

AARD, %

STD

10.266

0.0356

oo

3.3. Sensitivity Analyses over Wall’s Wettability

f

353

pr

Since it is known that wettability of the system greatly influences the capillary dominant

354 355 356

static contact angles. The dynamic contact angle is taken into consideration suing equation

357

(8).

358

Pr

e-

flow [66], in this section we investigate the saturation-Ca curves for tubes with different

al

Using equation (8), the advancing contact angle can be calculated by using static contact

359 360

(neutral), and 180 (oil-wet) were used for various Ca numbers investigated in this study. Note

361

that contact angle was defined based on the advancing fluid.

362

Jo u

rn

angle and Ca number as inputs. For this purposes, static contact angles of 0 (water-wet), 90

The calculated advancing contact angles are provided in Table 4. According to the

363

provided data, increasing Ca number in water-wet and neutral-wet mediums would result in

364

dynamic contact angles approaching 180 (i.e. oil wet medium). The only exception is the tube

365

with static contact angle of 180 in which dynamic contact angle will remain the same as the

366

static contact angle.

367 368

Table 4. Advancing contact angles calculated for displacements under different Ca numbers and three static contact angles

369 370

Journal Pre-proof

rn

al

Pr

e-

90 (Neutral)

oo

180 (Oil wet)

1.0E-4 1.0E-3 5.0E-3 8.0E-3 1.0E-2 2.0E-2 5.0E-2 8.0E-2 1.0E-1 1.0E+0 1.0E+1 1.0E-4 1.0E-3 5.0E-3 8.0E-3 1.0E-2 2.0E-2 5.0E-2 8.0E-2 1.0E-1 1.0E+0 1.0E+1 1.0E-4 1.0E-3 5.0E-3 8.0E-3 1.0E-2 2.0E-2 5.0E-2 8.0E-2 1.0E-1 1.0E+0 1.0E+1

Jo u

0 (Water wet)

Dynamic contact angle by Eq. (10) () 180 180 180 180 180 180 180 180 180 180 180 90.44 92.23 96.88 99.55 101.15 107.96 122.79 133.44 139.10 179.29 180.0 10.08 22.74 40.50 48.07 52.18 67.45 94.75 112.02 120.76 178.95 180.00

f

Ca number

pr

Static contact angle ()

Simulations were conducted by using the dynamic contact angles calculated at each Ca number and results are provided in the CDC curves in Fig. 5.

371 372 373

From the generated CDC curves, it is seen that there are three different regions for all of

374

the investigated static wettabilities. In the first low-Ca region, residual saturation is close to

375

zero at the time of breakthrough. In the second transition region, residual saturation increases

376

with increases in the Ca number. Eventually, in the third high-Ca region, residual saturation

377

Journal Pre-proof stabilizes in a certain value for all different wettabilities. The low-Ca and high-Ca regions are

378

the same for all the cases; the only difference between the different wettabilities occurs

379

within the transition zone.

380 381

the other wetting mediums. The transition zone starts almost at Ca=1E-3 for the oil-wetting

382

case, around Ca=1E-2 for the neutral-wetting case, and in Ca=2E-2 for the water-wetting

383

medium. As a result, the transition region for the oil-wetting medium begins at Ca numbers

384

f

The transition zone starts early for the oil-wetting medium and is smoother compared to

385

regions is almost the same for neutral wetting and water-wetting mediums. Finally, all

386

pr

oo

one log cycle smaller compared with other cases. The commencement of the transition

387

beginning of the transition zone is different for different cases, the end of the transition zone,

388

which happens in Ca=1E-1, seems to be same for all the cases. Saturation and pressure

389

profiles in the central plane have been provided to understand why transition zone is different

390

for each wettability.

391

Jo u

rn

al

Pr

e-

different wetting mediums converge to the residual saturation of around 50%. Although the

Figure 5. CDC curves generated for different fluid-wall static contact angles using VOFCSF scheme assuming dynamic contact angle conditions

392 393 394 395 396

Saturation map along the centerline of the tube is plotted for times close to the

397

breakthrough in Fig. 6. As it is seen, increasing the Ca number causes the displacing fluid to

398

flow faster along the center of the tube, thereby bypassing the fluid close to the wall. This

399

will generate a thick film of the displaced fluid that will remain attached to the wall and is

400

responsible for the increases in the residual saturation with increases in Ca number. However,

401

Journal Pre-proof 402

the saturation map as well. As it is seen for the neutral wettability in Ca=1E-2, the change in

403

contact angle from static value of 90 to the dynamic value of 101.15 (see Table 4) has led to

404

the separation in the displacing phase. This will in turn increase the residual saturation at the

405

time of breakthrough and is responsible for the early commencement of the transition zone in

406

the neutral-wet case compared to the water-wetting medium. At the high Ca number of 1.0,

407

all the cases exhibit similar saturation map since viscous force is the dominant force and also

408

all of them has reached the same dynamic contact-angle. As a result, the dynamic contact

409

Jo u

rn

al

Pr

e-

pr

angle might be only important for the transition zone.

oo

f

changing Ca number will also change the contact-angle of the fluid-wall which will influence

410 411 412 413 414 415 416 417 418 419 420 421 422 423

Journal Pre-proof 424

Static wettability

Ca

Saturation Map

1E-4 5E-3

Water- wet

oo

f

1E-2

pr

1.0

e-

1E-4 5E-3

Pr

Neutral-wet

Jo u

1E-4

rn

1.0

al

1E-2

5E-3

Oil-wet

1E-2 1.0

425

Figure 6. Saturation map (dark blue=100% water and red=100% oil) along the centerline for simulations under dynamic contact angle assumption

426 427 428

Journal Pre-proof To compare the velocity profiles for all the investigated cases, contours of velocity

429

magnitude are plotted in Fig. 7 for four different Ca numbers. The scale of the velocity

430

profiles is different for each case, but we just provided the qualitative legend to show the

431

distribution of the velocity magnitude along the central plane of the tube. The magnitudes of

432

the velocities are not important here; we are only interested in distribution of the velocity

433

magnitudes.

434

f

Although saturation profiles are similar for all the wettabilities in low-Ca region, it seems

435 436

happens right at the front for the water-wet system while for the neutral-wet system, the

437

pr

oo

that velocity profiles would be different for this region. In Ca=1E-4, the maximum velocity

438

maximum magnitude of the velocity for oil-wetting medium in Ca=1E-4 also happens near

439

the water-oil interface; however, it has more uniform velocity profile compared to the water-

440

wet medium.

441

al

Pr

e-

maximum velocity happens uniformly in the area extended from the inlet to the front. The

rn

Nonetheless, the velocity profile in the transition zone is almost the same for all the cases;

442 443

phase causes a non-uniform velocity map compared to other cases. Altering the Ca number to

444

the high-Ca zone, i.e. Ca=1.0, shrinks the width of the maximum velocity zone. Ultimately,

445

the velocity profile would be different for the low-Ca and transition regions, for mediums

446

with different static contact angles.

447

Jo u

the only exception is the neutral-wet case at Ca=1E-2 in which separation of the displacing

448 449 450 451

Journal Pre-proof

Static wettability

Ca

Saturation Map

1E-4 5E-3

Water- wet 1E-2

oo

f

1.0

pr

1E-4 5E-3

e-

Neutral-wet

Pr

1E-2

5E-3

Jo u

Oil-wet

rn

1E-4

al

1.0

1E-2 1.0

452

Figure 7. Velocity profile on the centerline for various Ca numbers and static contact angle of 180 (water wet)

453 454 455

Journal Pre-proof Finally, we generated “phase-diagram” of Lenormand for all three different wettabilities

456

to detect the flow-regimes involved in our simulations. In the phase-diagram, saturation of

457

the displacing fluid at breakthrough is plotted against the logarithm of the Ca number, for a

458

given viscosity ratio. The shape of the plotted curve will determine the types of the different

459

flow regimes [13].

460 461

cases exhibit “capillary fingering” and “viscous fingering” as two possible flow-regimes. To

462

f

The phase-diagrams are plotted in Fig. 8, and it is seen that all three different wettability

463

defined the new viscosity ratio “M” as the reciprocal of the viscosity ratio defined previously

464

pr

oo

make the diagrams consistent with Lenormand‟s definition of viscosity ratio [13], we have

465

fingering happens abruptly for the neutral-wet and water-wet cases compared to the oil-wet

466

case.

467

Jo u

rn

al

Pr

e-

in this paper. From the plots, it is also seen that shifting from capillary-fingering to viscous-

468

Jo u

rn

al

Pr

e-

pr

oo

f

Journal Pre-proof

469

Figure . Lenormand‟s phase-diagram for all three different wettabilities investigated in this study

470 471 472

Journal Pre-proof 473

3.4. Proposed Correlations As it was seen from the previous section, there are three different regions in the CDC

474

curves, low-Ca region, a transition zone, and a high-Ca region. In the low-Ca region, the

475

residual saturation at the breakthrough would be zero. Therefore, in this study, we

476

developed correlations to estimate the residual saturation in the transition and high-Ca

477

region, for three different static contact angle conditions.

478 479

oo

f

For neutral-wet and water-wet cases, residual saturation is functions of both Ca

480

trying variety of different functions:

481

)

(

))[

(

,

)

e-

(

]

482

(13)

483

Pr

(

pr

number and dynamic contact angle; therefore, we proposed the following function after

Since dynamic contact angle is itself functions of Ca number and static contact angle,

485

(

))[

rn

(

(

)

(

)

(

)

]

,

(14)

486 487

Jo u

)

al

see equation (8), the above equation may be rewritten as following: (

484

For the oil-wetting medium, since dynamic contact angle remains constant in the

488

CDC curves (see Table 4), following functionality was used by simply omitting the

489

cosine term from equation (13):

490

(

)

(

(

,

))

(15)

In the above equations, the terms a, b, and c are correlation coefficients which should be determined after applying the optimization algorithm.

491 492 493

We used nonlinear regression technique and employed Levenberg-Marquardt

494

technique as the optimization algorithm to find the correlation coefficients. Since

495

beginning of the transition zone is different for each case, the range of each developed

496

Journal Pre-proof correlation would be different as well. Within the Levenberg-Marquardt algorithm, sum

497

of square of errors was used as the cost function.

498 499

each correlation are provided in Table 5. Also, the match between the simulations and

500

regression data and the RE (Relative Error) for each Ca number are provided in figures 9-

501

11.

502

oo

f

The developed correlations for each medium and the range of the Ca numbers for

Static Wettability

pr

Correlation

Water wet

e-

Neutral wet

2×10-2 - 10 1×10-2 - 10 5×10-3 - 10

Pr

Oil wet

Ca Range

Jo u

rn

al

Table 5. The proposed correlations for different static wettability conditions

503 504 505 506 507 508 509 510 511

Journal Pre-proof 1.00

100

CFD Simulations (Water-wet) New correlation (Water-wet) RE

80

0.70

70

0.60

60

0.50

50

0.40

40

0.30

30

0.20

20

oo

f

Sor

0.80

90

Relative Error, %

0.90

0.10

1.0E-02

1.0E-01

1.0E+00

pr

0.00 1.0E-03

1.0E+01

10 0 1.0E+02

Ca

512 513

from the developed correlations for the oil-wet medium

514

Pr

e-

Figure 9. A comparison between the results of the CFD simulations and the results obtained

80 70

Jo u

Sor

0.70

90

0.60

60

0.50

50

0.40

40

0.30

30

0.20

20

0.10

10

0.00 1.0E-03

1.0E-02

1.0E-01

1.0E+00

1.0E+01

Relative Error, %

0.80

100

CFD Simulations (Neutral-wet) New correlation (Neutral-wet) RE

rn

0.90

al

1.00

0 1.0E+02

Ca

515

Figure 10. A comparison between the results of the CFD simulations and the results obtained

516

from the developed correlations for the Neutral-wet medium

517

Journal Pre-proof 100

CFD Simulations (Oil-wet) New correlation (Oil-wet) RE

0.90

80 70

0.60

60

0.50

50

0.40

40

0.30

30

0.20

20

f

0.70

oo

Sor

0.80

90

0.10

1.0E-02

1.0E-01

1.0E+00

pr

0.00 1.0E-03

1.0E+01

Relative Error, %

1.00

10 0 1.0E+02

e-

Ca

Pr

Figure 11. A comparison between the results of the CFD simulations and the results obtained

al

from the developed correlations for the Oil-wet medium

518 519 520 521 522

neutral-wet mediums with the maximum relative error of around 10%. For the neutral-

523

wetting and water-wetting mediums, the relative errors are scattered around the 0% error-

524

line while for the oil-wetting medium the RE vector is distributed around the 10% error-

525

line. In addition, the proposed correlations accurately capture the trends in the S or-Ca

526

curves and reach a plateau in high Ca numbers.

527

Jo u

rn

As it is seen from figures 9-11, excellent match is obtained for the water-wet and

The statistical parameters of AARD (%), correlation coefficient (R 2), and RMSE

528

(Root Mean Square Error) are calculated for the developed correlations and numerical

529

values are reported in Table 6. Accordingly, the maximum correlation coefficient is for

530

the water-wetting medium with the value of 0.9695. The neutral-wetting medium has a

531

high R2 of 0.94, and the oil-wetting medium still have a fair correlation coefficient of

532

Journal Pre-proof 0.8052. The maximum AARD among the investigated cases belongs to oil-wet case with

533

the value of around 12%. The minimum value of AARD=1.37% is for the water-wet case.

534

Also, the highest value of RMSE is 0.0461 for the oil-wetting medium, which is small

535

compared with the range [0,1] of the saturation parameter.

536 537

and RMSE indicate that the proposed correlations are accurate in estimating the residual

538

saturation of the displaced fluid within the transition and high Ca regions.

539

oo

f

Overall, excellent graphical match, higher values of R2 and lower values of AARD

Table 6. Statistical parameters for the developed correlations in this study R2

AARD, %

0.9695

1.37

Neutral wet

0.9406

3.58

Oil wet

0.8052

12.2

Minimum RE, % 0.122

Maximum RE, % 2.02

RMSE

0.69

9.61

0.0168

1.89

27.8

0.0461

0.0068

Pr

e-

pr

Static wettability Water wet

540

al

In closing, results of this study are applicable to the physical phenomena in which two-

rn

phase immiscible flow happens in thin capillary tubes. Examples of such fluid-fluid

Jo u

displacements can be observed in fuel cells and bioreactors [67], waterflooding EOR [68],

541 542 543 544

cementation in oil wells [69], flow in microchannel such as cracks and silts [70], capillary

545

trapping of CO2 in the CO2 sequestration process [71], flow in soils [72], and polymer

546

processing [63], provided that dimensionless number analyses of the flow lead to viscous or

547

capillary fingering flow regimes. The proposed correlations can be used in such processes

548

when there is an immediate need to estimate the residual saturation under operating

549

conditions, but there is not enough time to conduct experiments and/ or numerical

550

simulations.

551 552 553

Journal Pre-proof 554 555

In this study, we developed a numerical finite volume simulator using VOF-CSF strategy

556

to track the fluid-fluid interface and fluid-fluid-wall contact angle. Dynamic contact angle

557

was taken into considerations using an empirical equation. The developed model was verified

558

by using the experimental data adopted from the open literature, and sensitivity analyses were

559

conducted over Ca number and wall‟s wettability.

560

oo

f

4. Conclusion

561

distinct regions: low-Ca, transition, and high-Ca regions. The transition zone starts early for

562

the oil-wet case, compared to water-wet and neutral-wet cases. The saturation maps showed a

563

e-

pr

Results of this study unveiled that CDC curves of the residual saturation have three

Pr

separation of phases for the neutral-wet system in the transition zone. Furthermore, the

564 565

the same. Finally, three different correlations were developed to estimate the residual oil

566

al

velocity map was different for the low-Ca region even though saturation profiles remained

rn

saturation in transition and high-Ca regions for three different static wettabilities. The

Jo u

developed correlations showed high accuracy with the maximum AARD (%) of around 12%

567 568

among all the cases. Given the high accuracy of the developed correlations, they can be used

569

to estimate the residual oil saturation in the transition and high-Ca region when there is lack

570

of such data. Lenormand‟s phase-diagrams showed that the results of this study are applicable

571

to both capillary-fingering and viscous-fingering flow-regimes.

572 573

Acknowledgement: We would like to thank Dr. Howard A. Stone at Harvard University for his valuable comments on our paper.

574 575 576 577

Journal Pre-proof 578 579

Figure Captions:

580

correlations for residual saturation of the displaced fluid

581

Figure 2. The meshed capillary tube used for the numerical FV simulations

582

Figure 3. Average pressure at the inlet face, and the average oil saturation within the capillary tube versus time for coarse and adapted meshing geometries

583 584

f

Figure 1. A schematic of all the main steps utilized in this study to develop new

585 586

Figure 5. CDC curves generated for different fluid-wall static contact angles using VOF-

587

pr

oo

Figure 4. A comparison between the experimental and FV numerical results based on the saturation-Ca number curves

e-

CSF scheme and dynamic contact angle conditions

Pr

Figure 6. Saturation map along the centerline for simulations assuming dynamic contact angle conditions

al

Figure 7. Velocity profile on the centerline for various Ca numbers and static contact

rn

angle of 180 (water wet)

Jo u

Figure . Lenormand‟s phase-diagram for all three different wettabilities investigated in this study

Figure 9. A comparison between the results of the CFD simulations and the results obtained from the developed correlations for the oil-wet medium Figure 10. A comparison between the results of the CFD simulations and the results obtained from the developed correlations for the Neutral-wet medium Figure 11. A comparison between the results of the CFD simulations and the results obtained from the developed correlations for the Oil-wet medium

588 589 590 591 592 593 594 595 596 597 598 599 600 601

Journal Pre-proof 602

Nomenclatures

θ

603

Velocity vector

604

Contact angle

605

Pressure

606

Viscous stress tensor

607

Gravity acceleration

608

oo

f

̿

Fluid density

Body force due to surface tension

609

Color function

610

Time increment in the CFD simulation

611

e-

pr

F(σ)

Pr

Spatial increment in the CFD simulation Maximum of local Eigen values of a preconditioned

L So

n

rn

Nμ

Jo u

M

al

system

612 613 614

Viscosity ratio used in Lenormand‟s phase-diagram

615

Viscosity ratio

616

Average velocity of the injection front

617

Length of the capillary tube

618

Oil saturation

619

Surface tension

620

Curvature of the interface

621

Dirac delta function

622

The normal vector at the interface

623

The unit vector normal to the wall

624

The unit vector tangential to the wall

625

Journal Pre-proof Static advancing contact angle

626

Dynamic advancing contact angle

627

Sw

Water saturation

628

Sor,bt

Residual oil saturation at breakthrough

629 630 631

CFD: Computational Fluid Dynamics

632

VOF: Volume of Fluid

633

oo

f

Abbreviations:

CSF: Continuum Surface Force

pr

EOR: Enhanced Oil Recovery IFT: Interfacial Tension

Ca: Capillary number Re: Reynolds Number

Jo u

References:

rn

al

CDC: Capillary Desaturation Curve

Pr

e-

CFL: Courant-Friedrichs-Lewy

[1] E.J. Soares, R.L. Thompson, D.C. Niero, Physics of Fluids (1994-present) 27 (2015) 082105. [2] C. Han, Rheology, Springer, 1980, p. 121-128. [3] T. Kuppusamy, J. Sheng, J. Parker, R. Lenhard, Water Resources Research 23 (1987) 625. [4] J.F. Freitas, E.J. Soares, R.L. Thompson, International Journal of Multiphase Flow 37 (2011) 640. [5] E.J. Soares, P.R. Mendes, M.d.S. Carvalho, Journal of the Brazilian Society of Mechanical Sciences and Engineering 30 (2008) 160. [6] W. Olbricht, D. Kung, Physics of Fluids A: Fluid Dynamics (1989-1993) 4 (1992) 1347. [7] H. Geistlinger, I. Ataei-Dadavi, H.-J. Vogel, Transport in Porous Media 112 (2016) 207. [8] K. Rahimi, M. Adibifard, Petroleum Science and Technology 33 (2015) 79. [9] K. Rahimi, M. Adibifard, M. Hemmati, H. Shariat Panahi, S. Gerami, Asia‐Pacific Journal of Chemical Engineering 11 (2016) 98. [10] H. Goldsmith, S. Mason, Journal of Colloid Science 18 (1963) 237. [11] S. Iglauer, M. Fernø, P. Shearing, M. Blunt, Journal of colloid and interface science 375 (2012) 187. [12] R. Lenormand, C. Zarcone, A. Sarr, Journal of Fluid Mechanics 135 (1983) 337. [13] R. Lenormand, E. Touboul, C. Zarcone, Journal of fluid mechanics 189 (1988) 165. [14] S. Zendehboudi, N. Rezaei, I. Chatzis, Journal of Porous Media 15 (2012).

634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660

Journal Pre-proof 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682

engineering 29 (1981) 329. [33] H.K. Versteeg, W. Malalasekera, An introduction to computational fluid dynamics: the finite volume method. Pearson Education, 2007. [34] K.E. Wardle, T.R. Allen, R. Swaney, Separation Science and Technology 41 (2006) 2225. [35] A. Golpaygan, N. Ashgriz, International journal of computational fluid dynamics 22 (2008) 85. [36] Y. Ding, R. Anderson, L. Zhang, X. Bi, D.P. Wilkinson, International Journal of Multiphase Flow 52 (2013) 35. [37] K.E. Wardle, Separation Science and Technology 46 (2011) 2409. [38] K.E. Wardle, T.R. Allen, M.H. Anderson, R.E. Swaney, AIChE journal 54 (2008) 74. [39] A.B. Desamala, V. Vijayan, A. Dasari, A.K. Dasmahapatra, T.K. Mandal, Journal of Hydrodynamics, Ser. B 28 (2016) 658. [40] G.-y. LU, W. Jing, Z.-h. JIA, Journal of Hydrodynamics, Ser. B 19 (2007) 683. [41] J. Mencinger, B. Bizjan, B. Širok, International Journal of Multiphase Flow 77 (2015) 90. [42] M. Passandideh-Fard, E. Roohi, International Journal of Computational Fluid Dynamics 22 (2008) 97. [43] G.H. Schnerr, J. Sauer, Fourth international conference on multiphase flow, New Orleans, USA, 2001. [44] M.R. Davidson, M. Rudman, Numerical Heat Transfer: Part B: Fundamentals 41 (2002) 291. [45] C.W. Hirt, B.D. Nichols, Journal of computational physics 39 (1981) 201. [46] J. Jeong, D. Yang, International Journal for Numerical Methods in Fluids 26 (1998) 1127. [47] M. Dianat, M. Skarysz, A. Garmory, International Journal of Multiphase Flow 91 (2017) 19. [48] B. Van Wachem, A.-E. Almstedt, Chemical Engineering Journal 96 (2003) 81. [49] J. Brackbill, D.B. Kothe, C. Zemach, Journal of computational physics 100 (1992) 335. [50] A.A. Research, ANSYS, Inc., Release 16.0. [51] B.A. Nichita, I. Zun, J.R. Thome, 7th International Conference on Multiphase Flow, ICMF 2010, Tampa, FL, May 30–June 4, 2010, 2010. [52] R. Courant, K. Friedrichs, H. Lewy, Mathematische annalen 100 (1928) 32. [53] F. Fairbrother, A.E. Stubbs, Journal of the Chemical Society (Resumed) (1935) 527. [54] L.E. Brownell, D.L. Katz, Chemical Engineering Progress 43 (1947) 537.

683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711

Jo u

rn

al

Pr

e-

pr

oo

f

[15] I. Chatzis, N.R. Morrow, H.T. Lim, Society of Petroleum Engineers Journal 23 (1983) 311. [16] H. Westborg, O. Hassager, Journal of colloid and interface science 133 (1989) 135. [17] Z.-h. LIU, Y.-p. GAO, Journal of Hydrodynamics, Ser. B 19 (2007) 630. [18] J.Q. Feng, International Journal of Multiphase Flow 35 (2009) 738. [19] M. Adibifard, Physics of Fluids 30 (2018) 082107. [20] A. Nabizadeh, H. Hassanzadeh, M. Sharifi, M.K. Moraveji, Journal of Molecular Liquids (2019) 111457. [21] D.H. Rothman, Journal of Geophysical Research: Solid Earth 95 (1990) 8663. [22] A. Kovscek, H. Wong, C. Radke, AIChE Journal 39 (1993) 1072. [23] A.B. Dixit, J. Buckley, S.R. McDougall, K.S. Sorbie, Transport in Porous Media 40 (2000) 27. [24] M.J. Blunt, Journal of Petroleum Science and Engineering 20 (1998) 117. [25] A. Eslami, S. Taghavi, Journal of Non-Newtonian Fluid Mechanics 243 (2017) 79. [26] Y. Guo, X.B. Yu, International Journal of Multiphase Flow 91 (2017) 89. [27] A. de Lima Cunha, S. Neto, A.G.B. de Lima, E.S. Barbosa, Defect and Diffusion Forum, Trans Tech Publ, 2011, p. 935-940. [28] A. Nabizadeh, M. Adibifard, H. Hassanzadeh, J. Fahimpour, M.K. Moraveji, Journal of Molecular Liquids 273 (2019) 248. [29] T. Narasimhan, P. Witherspoon, Water Resources Research 12 (1976) 57. [30] E. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yusof, Journal of computational physicsJournal of computational physics 161 (2000) 35. [31] J. Marchal, M. Crochet, Journal of Non-Newtonian Fluid Mechanics 26 (1987) 77. [32] T.J. Hughes, W.K. Liu, T.K. Zimmermann, Computer methods in applied mechanics

Journal Pre-proof

Pr

e-

pr

oo

f

[55] E. Ojeda, F. Preston, J. Calhoun, Prod. Monthly 18 (1953) 28. [56] P.G. Saffman, G. Taylor, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, The Royal Society, 1958, p. 312-329. [57] W. Foster, Journal of Petroleum Technology 25 (1973) 205. [58] R. Ehrlich, H. Hasiba, P. Raimondi, Journal of Petroleum Technology 26 (1974) 1. [59] R.L. Reed, R.N. Healy, Improved oil recovery by surfactant and polymer flooding, Elsevier, 1977, p. 383-437. [60] J. Garnes, A. Mathisen, A. Scheie, A. Skauge, SPE/DOE Enhanced Oil Recovery Symposium, Society of Petroleum Engineers, 1990. [61] H. Guo, M. Dou, W. Hanqing, F. Wang, G. Yuanyuan, Z. Yu, W. Yansheng, Y. Li, SPE Kuwait Oil and Gas Show and Conference, Society of Petroleum Engineers, 2015. [62] R.A. Fulcher Jr, T. Ertekin, C. Stahl, Journal of petroleum technology 37 (1985) 249. [63] E.J. Soares, R.L. Thompson, D.C. Niero, Physics of Fluids 27 (2015) 082105. [64] H.K. Esfeh, A. Azarafza, M. Hamid, RSC Advances 7 (2017) 32893. [65] M.J. Berger, P. Colella, Journal of computational physics 82 (1989) 64. [66] S. Zendehboudi, N. Rezaei, I. Chatzis, Energy & Fuels 25 (2011) 4452. [67] P. Rapolu, S.Y. Son, ASME 2007 International Mechanical Engineering Congress and Exposition, American Society of Mechanical Engineers Digital Collection, 2009, p. 1039-1045. [68] Z. Kargozarfard, M. Riazi, S. Ayatollahi, Petroleum Science 16 (2019) 105. [69] E. Soares, M. Carvalho, P. Mendes, Journal of Fluids engineering 127 (2005) 24. [70] S. Ghiaasiaan, S. Abdel-Khalik, Advances in heat transfer, Elsevier, 2001, p. 145-254. [71] C.H. Pentland, R. El‐Maghraby, S. Iglauer, M.J. Blunt, Geophysical Research Letters 38 (2011). [72] E.W. Washburn, Physical review 17 (1921) 273.

Declaration of interests

rn

al

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

   

Jo u

Highlights

712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740

A numerical model of a micro-tube was generated using VOF-CSF strategy for fluid-

741

fluid and fluid-fluid-surface tracking

742

Adaptive Mesh Refinement (AMR) was used for Mesh sensitivity analysis and

743

experimental data used for physical validation of the model

744

CDC (Capillary Desaturation Curve) curves were generated considering dynamic

745

contact angle

746

Three different correlations were proposed for residual saturation of the displaced

747

fluid with maximum average relative error of 12%

748

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Figure 11