1-s2 0-S037702572200074X-main

Prévia do material em texto

Journal of Non-Newtonian Fluid Mechanics 305 (2022) 104834
A
0
F
d
P
D
A
K
F
Y
C
f
1
i
a
t
a
t
s
v
f
M
p
i
f
i
u
i
f
h
R
Contents lists available at ScienceDirect
Journal of Non-Newtonian Fluid Mechanics
journal homepage: www.elsevier.com/locate/jnnfm
low instabilities in fluid displacement through enlarged regions in annular
ucts
.R. Varges, E.C. Rodrigues, L.C. Moraes, P.R. de Souza Mendes, M.F. Naccache ∗
epartment of Mechanical Engineering, Pontifícia Universidade Católica do Rio de Janeiro, Rua Marquês de São Vicente 225, Rio de Janeiro, RJ
22453-900, Brazil
R T I C L E I N F O
eywords:
low displacement
ield stress fluids
omputational fluid dynamics
low instabilities
A B S T R A C T
This work presents numerical results of flow displacement and the instabilities that appear when a viscoplastic
fluid is displaced by a Newtonian one in a vertical annular duct with an enlarged region. Hydrodynamic
instabilities in displacement flows may be due to viscosity and density differences, inertia effects and geometry
irregularities like eccentricity and variable cross sections. If one of the fluids is non-Newtonian, the situation
can be more challenging due to the dependence of the rheology on the flow kinematics. The geometry analyzed
is found in the well cementation process in the oil industry, where the drilling fluid (viscoplastic behavior) is
displaced by a spacer fluid (Newtonian). The drilling fluid wash out must be perfect in order to guarantee the
success of the cementation operation, thus preventing the collapse of the oil well. Situations with excessive
wellbore enlargement can occur, hampering the removal of the drilling fluid and the establishment of a perfect
zonal isolation. The numerical solution of the problem is obtained using the finite volume technique to solve
the governing conservation equations of mass and momentum. The volume of fluid method is employed to deal
with the multiphase problem. The aim of this study is to look at how the governing parameters affect the flow
displacement pattern and efficiency inside the enlarged region, with the goal to optimize the displacement
process. The results focus on real case scenario for unstable situations, and present new findings that helps
understanding the role of viscous, inertia and buoyancy effects.
. Introduction
The displacement flow of a fluid by another has been widely studied
n the literature due to its relevance to several natural phenomena
nd industrial processes. Depending on the flow and fluid parame-
ers, hydrodynamic instabilities may appear. When a fluid displaces
more viscous one, viscous fingering occurs, a phenomenon called
he Saffman–Taylor instability. Density differences may also play a
ignificant role in the onset of instabilities, either for horizontal or
ertical displacements. The complexity increases when non-Newtonian
luids are present, due to viscosity dependence on the flow kinematics.
oreover, elasticity and thixotropy can affect the flow displacement
attern. Besides fluid properties, flow rate and geometry can play an
mportant role, often leading to instabilities that can improve or impair
low displacement. Depending on the application, the hydrodynamic
nstabilities can be desirable or not, but in any case it is essential to
nderstand which parameters are responsible for them, and how they
nfluence their development.
Displacement flows through abruptly varying cross sections are
ound in several industrial and natural processes, such as flows through
∗ Corresponding author.
E-mail address: naccache@puc-rio.br (M.F. Naccache).
URL: http://greo.mec.puc-rio.br (M.F. Naccache).
porous media, blood flow circulation, extrusion, and cementation of
oil wells. The analysis presented here focus on the last one, in cir-
cumstances where an enlarged region is present. In this and other
situations, the fluids can present a non-Newtonian behavior, which can
lead to some special features that result in a more complex analysis. For
example, when a viscoplastic fluid is displaced, it will remain stagnant
in the regions where the stress do not surpass the yield stress, resulting
in a non-effective displacement and changing significantly the pressure
drop throughout the flow. The most probable outcome will be a non
perfect final product or process. These issues are critical in the oil
industry, not only during cementing and drilling operations but also
in the displacement through reservoirs during production [1–5]. In
several other areas, flow displacement problems appear in flows during
extrusion processes. Another example of application, in the area of
biomedicine, is studied by Roberts et al. [6], who analyzed the process
of using flow of a viscoplastic aqueous foam to deliver surfactant to
a varicose vein, and showed that increasing the yield stress leads to a
more effective blood displacement inside the veins.
vailable online 25 May 2022
377-0257/© 2022 Elsevier B.V. All rights reserved.
ttps://doi.org/10.1016/j.jnnfm.2022.104834
eceived 29 October 2021; Received in revised form 17 March 2022; Accepted 17
 May 2022
http://www.elsevier.com/locate/jnnfm
http://www.elsevier.com/locate/jnnfm
mailto:naccache@puc-rio.br
http://greo.mec.puc-rio.br
https://doi.org/10.1016/j.jnnfm.2022.104834
https://doi.org/10.1016/j.jnnfm.2022.104834
http://crossmark.crossref.org/dialog/?doi=10.1016/j.jnnfm.2022.104834&domain=pdf
Journal of Non-Newtonian Fluid Mechanics 305 (2022) 104834P.R. Varges et al.
Fig. 1. Schematic of the abrupt annular expansion–contraction geometry.
The effects of flow rate and fluid rheology through regular and irreg-
ular geometries have been analyzed for single fluids both numerically
and experimentally [7–15]. These researchers studied the role of flow
rate and yield stress in the appearance and development of unyielded
regions that directly affect the flow pattern and pressure drop. They
observed the presence of unyielded stagnant zones near corner regions
and/or unyielded plug zones close to the center of the duct.
Regarding the displacement of one fluid by another, not only the
fluids rheology, flow rate and geometry are relevant, but also and more
importantly the viscosity and density differences. Flow displacement in
annular spaces has been analyzed in the literature focusing mainly on
cementing operations of oil wells [15–23], where it has been examined
the effects of viscosity and density ratios, flow rate and annulus eccen-
tricity on the displacement efficiency. Geometry irregularities are also
addressed, due to its importance for the oil industry. The cementing
operation in a vertical washout section of a well (annular duct with an
expansion followed by a contraction) has also been analyzed [1–3].
The flow of miscible viscoplastic fluids in eccentric annuli was
analyzed numerically and experimentally by Deawwanich et al. [20].
Flow images showed the effects of eccentricity, cylinder rotation and
rheology on the displacement. Etrati et al. [2] presented a numerical
study of the displacement of a Bingham fluid (modeling the drilling
fluid) by a Newtonian one (the spacer fluid) to investigate the effect
of the washout length, flow rate, density difference and viscosity ratio
2
on the displacement efficiency. The results are limited to positive
density differences (heavier fluid displacing a lighter one) and higher
Reynolds numbers (but still within the laminar range), showing that
density differences play a major role in the process, with larger density
differences improving displacement efficiency. The authors also ob-
served unstable situations as the Reynolds number increases, but they
verified that the instabilities may improve or impair the displacement.
Some unyielded fluid is left behind, which is very damaging in the
case studied. Naccache et al. [1] also performed a numerical study,
using different combinations of pairs of fluids. The results showed the
transition from viscous-dominated (low Re), where a fore–aft symmetry
is observed inside the washout zone, to inertia-dominated(higher Re)
regimes, where the symmetry breaks and the interface between fluids
is displaced to the wall near the contraction zone. Flow rate, rheology
and geometry effects on the displacement efficiency were investigated.
Espinoza et al. [3] presented experimental and numerical results in a
concentric annuli washout zone in a low scale geometry. All the cases
analyzed resulted in stable situations, and it is observed that increasing
the Reynolds number, washout depth and length, the displacement
improves. However, the viscosity ratio did not seem to be relevant.
Renteria et al. [24] conducted experiments and 2D numerical sim-
ulations to study eccentricity and inclination effects on displacement
flows of yield stress materials in near-horizontal irregular annuli. It was
observed a poor displacement efficiency when considering fluids with
high yield stresses flowing in eccentric annuli at 0 or 80◦. Moreover,
it was noticed that when the displacement in the annulus is already
poor, enlarged regions worsen the displacement efficiency. On the
other hand, when the regular annulus presents high displacement effi-
ciency, the enlarged region may improved the displacement. Skadsem
et al. [25] also investigated experimentally and numerically the effects
of eccentricity and inclination in an irregular annulus characterized
by an enlarged region. The variation in axial flow velocity across the
cross-section of the annulus was observed to be the most prominent
consequence of eccentricity. The difference in flow velocity caused an
elongated interface in the regular part of the annulus, as well as a
poor displacement in the horizontal annulus enlargement. Furthermore,
they observed that as inclination decreases, the displacement efficiency
increases within concentric and eccentric annuli and within the regular
part of the annulus in front of the enlarged region.
Another experimental study of interest was performed by Akbari
and Taghavi [26], who analyzed the process of placement of a heavy
fluid to remove a lighter fluid in an inclined closed-ended pipe. This
process simulates a method used for plug and abandonment opera-
tions of oil wells. The two fluids were miscible and Newtonian, with
different densities. The results showed instabilities and mixing in the
front interface. It was shown that density and viscosity differences
increase the placement efficiency. Instabilities in flow displacement
were also observed by Tehrani et al. [27]. The authors analyzed the
flow through upward inclined annular eccentric cylinders, and showed
the channeling of the displacing fluid in the wide side of the geom-
etry, worsening the displacement process. Density difference (denser
displacing fluid) promotes azimuthal flow of the denser fluid towards
the narrow side, which, depending on the flow rate (or on the ratio
between the buoyancy to Reynolds number), may lead to instabilities
that improve the displacement.
Fig. 2. Volume fraction for 𝑡∗ = 1, 𝐷∗ = 1.5, 𝑑′∗ = 0.5, 𝑙∗ = 0.5, 𝐿∗ = 6.4, 𝑅𝑒 = 3.92, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.72, 𝑛 = 0.78, 𝛥𝜌∗ = 0, and 𝐺𝑎 = 0. Flow from the left (bottom) to right (top).
(a) Mesh 1 with 295,554 nodes; (b) Mesh 2 with 348,824 nodes; (c) Mesh 3 with 401,112 nodes. Fluid 2 is represented by black color and fluid 1 by gray color.
Journal of Non-Newtonian Fluid Mechanics 305 (2022) 104834P.R. Varges et al.
Fig. 3. Velocity profile at the center of the enlarged region for 𝑡∗ = 1, 𝐷∗ = 1.5,
𝑑′∗ = 0.5, 𝑙∗ = 0.5, 𝐿∗ = 6.4, 𝑅𝑒 = 3.92, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.72, 𝑛 = 0.78, 𝛥𝜌∗ = 0, and
𝐺𝑎 = 0 for the meshes with 295,554, 348,824, and 401,112 nodes.
Renteria et al. [28] and Sarmadi et al. [29] explored borehole
irregularities in primary cementing, specially in horizontal wells in-
volving yield stress materials. The first [28] focus on displacement
flows in eccentric and elliptic geometries, while Sarmadi et al. [29]
focus on circular geometries with enlargements of different depths and
wavelengths. It was observed that borehole irregularities can reduce the
displacement efficiency. Additionally, Jung and Frigaard [30] recently
investigated common cementing practices in horizontal wells to im-
prove the displacement efficiency. Excess of cement volume, centralizer
usage and staged cementing operations were numerically evaluated and
fully discussed.
Truzzolillo and Cipelletti [31] presented a review of the general
mechanisms driving hydrodynamic instabilities in miscible fluids, ad-
dressing the differences between miscible and immiscible interfaces.
Atmakis and Kenig [32] presented numerical results of the Kelvin–
Helmholtz (KH) instability in a 2D model geometry of two liquid
layers, using the Volume of Fluid method (VoF) and level-set methods.
The authors observed the stabilization effect of the surface tension
over the unstable finger-like structures. Lee and Kim [33] studied
numerically the Kelvin–Helmholtz instability of multi-component fluids
using a phase-field model, investigating the effects of surface tension,
density ratio, and magnitude of velocity difference. It was shown that
increasing the surface tension or the density ratio reduces the growth
of the KH instability. It was also observed that as the initial horizontal
velocity difference gets larger, the interface rolls up more. Furthermore,
Dihn et al. [34] observed that the presence of viscosity slows the rate
at which the interface grows.
In the present work, we extend the work of Espinoza et al. [3], using
real data of cementing operations for the geometry of the enlarged
3
Table 1
Mesh sizes, displacement efficiency 𝜙 and its difference with the most refined mesh
result 𝜖.
Mesh Number of nodes 𝜙 % 𝜖 %
Mesh 1 295,554 85.80 1.27
Mesh 2 348,824 85.22 0.59
Mesh 3 401,112 84.72 0
zone, fluids rheology and flow rates. It is worth mentioning that the
experimental and numerical results presented in [3] are all related to
stable situations and lower flow rates than the real cases. Therefore,
the present work aims to contribute to a better understanding of the
role of viscosity and density differences, geometry and inertia in the
upward flow displacement through expansion–contraction flows, and
how these parameters affect the displacement efficiency. Although
the main motivation of this study is to investigate the displacement
flow to optimize washout cementing operations in the oil industry,
the results can be used to avoid/induce instabilities to increase the
displacement efficiency of other similar processes, since all results are
analyzed via dimensionless parameters. We present numerical results
of a Newtonian fluid displacing a yield stress material modeled by a
regularized Herschel–Bulkley viscosity function. The effects of some
rheological, flow and geometric parameters on flow pattern/volume
fraction distribution and displacement efficiency are presented and dis-
cussed. Understanding the role of the governing parameters in the onset
and development of the hydrodynamic instabilities should provide
means to the optimization of the flow displacement process.
2. Modeling
Scheme of the flow geometry is presented in Fig. 1. The tube is
positioned vertically, and the upstream and downstream length and
diameter of the outer tube are equal to 𝑙 and 𝑑, respectively. The
diameter of the inner tube is 𝑑′, and the length and diameter of the
enlarged region are 𝐿 and 𝐷, respectively. Initially all the annular space
is filled with fluid 1. At the beginning of the simulation, fluid 2 enters
the tube from the bottom. Fluids are considered incompressible, fluid
2 is Newtonian and fluid 1 is a yield stress material.
Flow is isothermal and axisymmetric, since only cases of concentric
inner and outer tubes are investigated.
The fluids are immiscible, and the volume of fluid (VoF) method
is used to track the interface between the fluids [35]. This method
calculates the volume fraction of each phase 𝑗, 𝛼𝑗 = ∀𝑗∕∀ (𝑗 = 1,
2, . . . , n), where 𝑛 is the number of phases and ∀𝑗 is the volume
occupied by phase𝑗 in the cell control volume ∀. In the problem
analyzed, 𝑛 = 2 so that the volume fractions are obtained solving the
mass conservation equation for one phase plus a restriction equation,
Eqs. (1), (2), respectively given by:
𝜕𝛼1
𝜕𝑡
+ 𝐯 ⋅ ▽𝛼1 = 0 (1)
𝛼 + 𝛼 = 1 (2)
1 2
Fig. 4. Volume fraction at 𝑡∗ = 1 for 𝐿∗ = 6.4, 𝐷∗ = 1.5, 𝛥𝜌∗ = 0, 𝐺𝑎 = 0, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78 and: (a) 𝑅𝑒 = 0.04; (b) 𝑅𝑒 = 3.88; (c) 𝑅𝑒 = 93.1. Fluid 2 is represented by
black color and fluid 1 by gray color.
Journal of Non-Newtonian Fluid Mechanics 305 (2022) 104834P.R. Varges et al.
𝜂
w
z
a
c
a
b
𝛷
t
▽
𝜌
I
a
g
a
c
c
h
v
2
l
𝐷
𝜏
d
Fig. 5. Displacement efficiency through time for 𝐿∗ = 6.4, 𝐷∗ = 1.5, 𝛥𝜌∗ = 0, 𝐺𝑎 = 0,
∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78 and: (a) 𝑅𝑒 = 0.04; (b) 𝑅𝑒 = 3.88; (c) 𝑅𝑒 = 93.1.
here v is the velocity vector. The volume fraction is equal to one or
ero when a control volume is entirely filled with one of the phases,
nd it is equal to a value between one and zero if the control volume
ontains both fluids. The velocity and pressure at each control volume
re equal for both phases, while the properties are given by the average
etween the properties (𝛷) of each fluid:
= 𝛼1𝛷1 + (1 − 𝛼1)𝛷2 (3)
The mass and momentum conservation equations, and the constitu-
ive equation are given by:
⋅ 𝐯 = 0 (4)
[ 𝜕𝐯
𝜕𝑡
+ (𝐯 ⋅ ▽)𝐯
]
= −▽𝑝 +▽ ⋅ 𝝉 + 𝜌𝐠 + 𝐅𝐬 (5)
𝝉 = 2𝜂(�̇�)𝐃 (6)
In the above equations, 𝜌 is the density, 𝑝 is the pressure, 𝐠 is
the gravity vector, 𝐃 ≡ 1∕2
[
▽𝐯 +▽𝐯𝑇
]
is the rate-of-strain tensor,
�̇� ≡
√
2tr𝐃2 is the magnitude of 𝐃, and 𝐅𝑠 = 𝜎𝛿1▽𝛼1 is the surface
force per unit volume [36], where 𝜎 is the interfacial tension between
the two fluids and 𝛿1 is the curvature of the interface of fluid 1, defined
as the divergence of the unit normal vector. 𝝉 is the deviatoric stress
tensor, given by the Generalized Newtonian Fluid constitutive equation
(Eq. (6)). The displacing fluid, fluid 2, is Newtonian, so the viscosity is
constant and equal to 𝜇2. The displaced fluid (fluid 1) is a yield stress
material, modeled by a regularized Herschel–Bulkley equation [37],
given by:
⎧
⎪
⎨
⎪
⎩
𝜂(�̇�) = 𝜏𝑦
�̇� + 𝑘�̇�𝑛−1, if �̇� > �̇�cr
𝜂(�̇�) = 𝜏𝑦
�̇�cr
[
2 − �̇�
�̇�cr
]
+ 𝑘�̇�𝑛−1cr
[
(2 − 𝑛) + (𝑛 − 1) �̇�
�̇�cr
]
, otherwise.
(7)
n the above equation, 𝜏𝑦, 𝑘 and 𝑛 are the yield stress, consistency index
nd power-law index of fluid 1. �̇�cr is the regularization parameter,
iven by the small rate-of-strain below which the viscosity is equal to
constant high value ≈ 500𝜂𝑐 . Therefore, �̇�cr ≈ 𝜏𝑦
500𝜂𝑐
, where 𝜂𝑐 is a
haracteristic viscosity. In the present work 𝜂 = 𝜂(�̇� ), where �̇� is a
4
𝑐 𝑐 𝑐 𝑙
haracteristic shear rate defined as �̇�𝑐 = 𝑉 ∕𝐷ℎ, 𝐷ℎ = 𝐷 − 𝑑′ being the
ydraulic diameter of the central (enlarged) tube, and 𝑉 is the mean
elocity at this enlarged region.
The initial and boundary conditions are:
• At the instant of time t = 0, the tube is filled with fluid 1 (𝛼1 = 1
and 𝛼2 = 0), and the velocity vector is zero (𝐯 = 𝟎).
• At the tube entrance, the Newtonian fluid (fluid 2) is injected with
a developed (parabolic) flow profile (𝐯 = 𝑢𝐞𝑥), with an average
velocity 𝑣𝑖𝑛.
• At the outlet, the flow is assumed fully developed.
• At the walls, no slip and impermeability conditions are applied, so
that both axial and radial velocity components are equal to zero.
.1. Scaling
The flow under analysis is scaled based on the following dimension-
ess variables.
∗ = 𝐷
𝐷ℎ
𝑑′∗ = 𝑑′
𝐷ℎ
𝑑∗ = 𝑑
𝐷ℎ
𝑙∗ = 𝑙
𝐷ℎ
𝐿∗ = 𝐿
𝐷ℎ
𝑡∗ = 𝑡
𝑡𝑐
∗
𝑦 =
𝜏𝑦𝐷ℎ
𝜇2𝑉
𝛥𝜌∗ =
𝜌2 − 𝜌1
𝜌2
𝜂∗ =
𝜇2
𝜂𝑐
𝜎∗ = 𝜎
𝜇2𝑉
(8)
The dimensionless time 𝑡∗ is obtained using the characteristic time
(𝑡𝑐) that represents the instant of time when exactly one domain volume
of fluid 2 has crossed the inlet boundary. Besides the dimensionless
geometric quantities (𝐷∗, 𝑑′∗, 𝑑∗, 𝑙∗, 𝐿∗), the dimensionless time 𝑡∗, the
dimensionless yield stress 𝜏∗𝑦 , the dimensionless density difference 𝛥𝜌∗,
the viscosity ratio 𝜂∗, and the dimensionless interfacial tension 𝜎∗,
some relevant dimensionless parameters can be defined to address
the balance among inertia, viscous and buoyancy forces. Then, the
Reynolds and the Galilei numbers are respectively defined as:
𝑅𝑒 =
𝜌2𝑉 𝐷ℎ
𝜇2
(9)
𝐺𝑎 =
(𝜌2 − 𝜌1)𝑔𝐷2
ℎ
𝜇2𝑉
(10)
To evaluate quantitatively the flow displacement we define the
isplacement efficiency 𝜙 as
𝜙 =
∀2
∀𝑇
(11)
where ∀2 is the volume occupied by the displacing fluid (fluid 2) at an
instant of time and ∀𝑇 is the total volume of the geometry.
3. Numerical solution
The conservation equations of mass and momentum are solved using
Fluent software (Ansys Inc.). The finite volume method is used to
discretize the conservation equations, and the volume of fluid method
was employed to deal with the multiphase flow. The pressure–based
solver is used for the transient model, and the PISO algorithm is
employed to couple the pressure–velocity equations. The mass conser-
vation equation is discretized with the power law method, and the
upwind discretization method is employed for the momentum equation
[19].
A validation of the numerical procedure with experimental results
was performed in a previous work [3], using different geometries and
parameters. The results are presented in the Supplementary Material
(Fig. S1), where we also show comparisons of 2D and 3D simulations
(Fig. S2). In all simulations analyzed in the present work, the flow
is assumed to be axisymmetric, so only 2D (𝑟 − 𝑥) situations were
investigated. A structured non-uniform mesh was generated using the
open source Gmsh software. Mesh tests were performed for a base case,
whose values of the governing parameters are 𝐷∗ = 1.5, 𝑑′∗ = 0.5,
∗ ∗ ∗ ∗ ∗
= 0.5, 𝐿 = 6.4, 𝑅𝑒 = 3.92, 𝛥𝜌 = 0, 𝜂 = 0.08, 𝜏𝑦 = 12.7, 𝑛 = 0.78.
Journal of Non-Newtonian Fluid Mechanics 305 (2022) 104834P.R. Varges et al.
Fig. 6. Volume fraction at 𝑡∗ = 1 for 𝑅𝑒 = 1, 𝜂∗ = 0.02, 𝜏∗𝑦 = 47.7, 𝑛 = 0.78, 𝐷∗ = 1.5, 𝐿∗ = 6.4, and : (a) 𝛥𝜌∗ = −0.67 and 𝐺𝑎 = −2.4 × 105, (b) 𝛥𝜌∗ = −0.01 and 𝐺𝑎 = −3.5 × 103 (c)
𝛥𝜌∗ = 0 (isodense) and 𝐺𝑎 = 0, (d) 𝛥𝜌∗ = 0.13 and 𝐺𝑎 = 4.7 × 104. Fluid 2 is represented by black color and fluid 1 by gray color.
Fig. 7. Displacement efficiency through time for 𝑅𝑒 = 1, 𝜂∗ = 0.02, 𝜏∗𝑦 = 47.7, 𝑛 = 0.78,
𝐷∗ = 1.5, 𝐿∗ = 6.4, and : (a) 𝛥𝜌∗ = −0.67 and 𝐺𝑎 = −2.4 × 105, (b) 𝛥𝜌∗ = −0.01 and
𝐺𝑎 = −3.5 × 103 (c) 𝛥𝜌∗ = 0 (isodense) and 𝐺𝑎 = 0, (d) 𝛥𝜌∗ = 0.13 and 𝐺𝑎 = 4.7 × 104.
Three different meshes were analyzed and the results were compared.
Table 1 presents the meshes used and the displacement efficiency for
each case. The difference between the displacement efficiency at a
given mesh (𝜙𝑖) and the one obtained with the finest Mesh 3 (𝜙𝑀3),
defined as 𝜖 = (𝜙𝑖 − 𝜙𝑀3)∕𝜙𝑀3, is also presented. All meshes are
slightly concentrated close to the walls, to the entrance and to the cross
section transition regions, as it can be seen in detail in Fig. S3 of the
Supplementary Material. Fig. 2 shows the volume fraction field and
Fig. 3 shows the velocity profile at the center of the eroded region
for the three different meshes. Considering the differences obtained
and the simulation time, the mesh with 348,824 nodes was chosen
for this geometry. For this case, the minimum face area normalized by
the area of the eroded zone, equal to 𝐿(𝐷 − 𝑑), is equal to 5.6 × 10−5,
and the maximum is 2.8 × 10−3. Mesh tests were not performed for the
geometries with different 𝐿∗. Instead, we kept the average size of the
face areas at the same order of magnitude, for all geometries analyzed.
The time step used in all simulations was adaptive, keeping the
numerical Courant number 𝐶𝑜 = (𝛥𝑡∕𝛥𝑥)𝑣 less than 1, where 𝛥𝑡 and
5
𝛥𝑥 are the time step and the local longitudinal cell size, and 𝑣 is the
local velocity. In all cases, the time step wasaround 10−3 to 10−2
s throughout the simulation. The convergence criteria used at each
time step was to consider the scaled residual [37] for the mass and
momentum equations below 10−5. It is worth mentioning that the
simulations are very time consuming. A regular case pumping one total
volume of fluid 2 (𝑡∗ = 1) with the mesh chosen, took around 96 h
running in a PC with Intel core i7 processor (3.4 GHz, 16 GB) with
only 4 physical cores.
4. Results and discussion
The results of flow displacement were obtained for different com-
binations of flow parameters, fluids rheology and geometry. We focus
the study at the enlarged region, and analyze snapshots of the flow
displacement and results of the displacement efficiency through time,
to study the influence of the parameters on the stability and efficiency
of the flow displacement. Instabilities are expected to occur due to
buoyancy effects (when upper fluid is heavier than lower one), adverse
viscosity ratios (less viscous fluid displacing a more viscous one) and
inertia effects (above a critical Reynolds number). We begin showing
the results of an iso-dense base case. Then, we discuss the buoyancy
effects, including a density difference 𝛥𝜌∗. We continue investigating
the influence of Reynolds number 𝑅𝑒, viscosity ratio 𝜂∗, enlarged region
length 𝐿∗, and finally interfacial tension 𝜎∗.
In all cases, the dimensionless diameter of the enlarged region is
kept equal to 𝐷∗ = 1.5. The inner tube, and entrance and exit tube
dimensionless diameters 𝑑′∗ = 0.5 and 𝑑∗ = 0.6 were also held fixed.
The entrance and exit tubes length were kept equal to 𝑙∗ = 5(𝑑∗−𝑑′∗) =
0.46, chosen to avoid any disturbances inside the enlarged region due
to entrance and exit boundaries, while minimizing the computational
domain. Moreover, all the results were obtained with interfacial tension
equal to zero (𝐅𝐬 = 𝟎, Eq. (5)), except in the last subsection, when the
effect of interfacial tension is analyzed. The results were explored for a
range of parameters variation, which were based in real case scenario of
cementing operations representing a spacer fluid displacing the drilling
mud. The dimensionless density difference ranges from −0.7 to 0.13,
Reynolds number ranges from 0.04 to 100, the viscosity ratio varied
from 0.0008 to 0.08. The dimensionless enlarged region length range
was 1 to 14 and the dimensionless interfacial tension varies up to 0.2.
It is worth mentioning that all cases are in the laminar regime for both
fluids and in all regions (entrance and exit tubes and enlarged region).
4.1. Effect of density difference
We first present a result of an iso-dense base case (𝛥𝜌∗ = 0 and
𝐺𝑎 = 0), using the rheology of fluids employed in real cementing
Journal of Non-Newtonian Fluid Mechanics 305 (2022) 104834P.R. Varges et al.
Fig. 8. Volume fraction at 𝑡∗ = 1 for 𝐿∗ = 6.4, 𝐷∗ = 1.5, 𝛥𝜌∗ = −0.01, 𝐺𝑎 = −955.4, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78, and: (a) 𝑅𝑒 = 0.04, (b) 𝑅𝑒 = 3.88, and (c) 𝑅𝑒 = 93.1. Fluid 2 is
represented by black color and fluid 1 by gray color.
Fig. 9. Displacement efficiency through time for 𝐿∗ = 6.4, 𝐷∗ = 1.5, 𝛥𝜌∗ = −0.01,
𝐺𝑎 = −955.4, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78, and: (a) 𝑅𝑒 = 0.04, (b) 𝑅𝑒 = 3.88, and (c)
𝑅𝑒 = 93.1. The blue square points refer to 𝑅𝑒 = 3.88 at 𝑡∗ = 2 and 3.
operations of the oil well industry. For these cases, the enlargement
length is 𝐿∗ = 6.4. Fluid 1 represents the drilling fluid, modeled as
a yield stress material with the viscosity given by the regularized HB
equation (7), with 𝜏∗𝑦 = 47.7 and 𝑛 = 0.78, and the displacing fluid (fluid
2) is Newtonian, representing the spacer fluid that flows prior to the
cement paste. The viscosity ratio is kept equal to 𝜂∗ = 0.02. The results
of the volume fraction field for three different Reynolds numbers at 𝑡∗ =
1 are depicted in Fig. 4, while the displacement efficiency evolution
is shown in Fig. 5. Both figures show a very efficient displacement in
all three cases. However, the evolution of the displacement is different
among them. In the first two cases, the Reynolds number is low and
so viscous forces dominate over inertia. As a result, and due to the
low viscosity ratio, the flow develops with a large finger of fluid 2,
which occupies almost all the enlarged region in the end, but no
instabilities are verified. It can be observed that only a small layer of
the displaced fluid remains close to the enlarged region walls and close
to the entrance and exit corners. In these regions, the stress does not
overcome the yield stress, and fluid 1 remains static. Moreover, the
displacement is almost fore–aft symmetric for this range of 𝑅𝑒, with
a slight deviation at the entrance due to viscous effects. Comparing
both low 𝑅𝑒 cases, we can also observe that as 𝑅𝑒 increases, the level
of stress increases and the layer of fluid close to the walls is slightly
6
reduced. As Reynolds increases further (Fig. 4c) the inertia effect can
be noted. First, the symmetry is reduced, and the entrance corner is
filled with the displacing fluid, while an amount of fluid 1 remains
close to the exit corner. Moreover, we can observe small spots of the
displaced fluid all over the enlarged region. This is in fact a result of
the evolution of the flow displacement. Despite the fact that the dis-
placement efficiency is quite similar for all cases, the flow displacement
evolution is completely different. This is shown in Figs. S4 and S5 of the
Supplementary Material. When inertia plays a role, the flow destabilize,
and recirculating vortices appear and develops throughout the enlarged
region, as fluid 2 flows. The unstable flow promotes mixing between
fluids and results in an efficient displacement, although some very small
spots of the displaced fluid seem to remain in the middle of the enlarged
region. It is worth mentioning that the presence of fluid 1, mixed with
fluid 2 throughout the domain or close to the walls, can impair the
cementing operation so that both situations should be avoided. Fig. 5
also shows that the displacement efficiency increases linearly with time,
and reaches an asymptote before 𝑡∗ = 1 for the stable flows, meaning
that the remaining fluid close to the walls and corners will remain
static, even for larger values of pumped volume. For the higher 𝑅𝑒 case,
this plateau is not reached because of the instabilities that promote
fluids mixing.
Density difference plays a significant role on the onset of instability.
When fluid 2 is lighter than fluid 1 the density difference is negative
and a destabilizing mechanism occur. Fig. 6 shows the volume fraction
for the same geometry of the previous cases, for 𝑅𝑒 = 1, 𝜂∗ = 0.02,
𝜏∗𝑦 = 47.7, 𝑛 = 0.78, and (a) 𝛥𝜌∗ = −0.67 (𝐺𝑎 = −2.4 × 105), (b)
𝛥𝜌∗ = −0.01 (𝐺𝑎 = −3.5 × 103), (c) 𝛥𝜌∗ = 0 (𝐺𝑎 = 0), and (d)
𝛥𝜌∗ = 0.13 (𝐺𝑎 = 4.7 × 104). All cases have low viscosity ratio and
low 𝑅𝑒, so inertia is negligible. Moreover, except for the iso-dense
case, 𝐺𝑎 is high in modulus, which indicates that buoyancy forces
dominates over viscous ones. Figs. 6a and b show the cases where the
displaced fluid is heavier than the displacing one. These cases show a
completely different pattern from the iso-dense cases showed above.
The displacement is highly unstable and its efficiency is very poor,
which can also be observed in Fig. 7. In both cases, buoyancy forces
are responsible for the onset of instabilities, since inertia is negligible.
The heavier fluid 1 tends to remain at the lower region of the enlarged
zone. As fluid 2 enters the enlarged zone, it generates a recirculation
of fluid 1, which on the other hand carries pieces of fluid 2 throughout
the enlarged zone. The same snapshots are shown in color at the
Supplementary Material (Fig. S6). For the highest density difference
(𝛥𝜌∗ = −0.67), we note a three layer pattern: at the lower zone of the
enlarged region, close to the entrance, fluid 1 remains almost static,
at the upper zone, close to the exit, a layer of fluid2 (lighter) forms,
and between then there is a recirculation zone filled mostly with fluid
1. This can be better observed with the aid of the velocity magnitude
field, shown in the Supplementary Material (Fig. S7). As the density
difference decreases (𝛥𝜌∗ = −0.01) the mixing region dominates all over
the enlarged zone, but the velocities are lower because the driving force
Journal of Non-Newtonian Fluid Mechanics 305 (2022) 104834P.R. Varges et al.
Fig. 10. Time evolution of the volume fraction field for 𝑅𝑒 = 3.88, 𝐿∗ = 6.4, 𝛥𝜌∗ = −0.01, 𝐺𝑎 = −955.4, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78, and: (a) 𝑡∗ = 0.05; (b) 𝑡∗ = 0.10; (c) 𝑡∗ = 0.15;
(d) 𝑡∗ = 0.2; (e) 𝑡∗ = 0.25; (f) 𝑡∗ = 0.5; (g) 𝑡∗ = 0.75; (h) 𝑡∗ = 1. Fluid 2 is represented by black color and fluid 1 by gray color.
Fig. 11. Volume fraction for: 𝑅𝑒 = 3.88, 𝐿∗ = 6.4, 𝛥𝜌∗ = −0.01, 𝐺𝑎 = −955.4, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78, 𝑡∗ = 1, and: (a) 𝑡∗ = 2 and (b) 𝑡∗ = 3. Fluid 2 is represented by black
color and fluid 1 by gray color.
Fig. 12. (a) Velocity intensity and (b) strain rate fields at 𝑡∗ = 1 for: 𝑅𝑒 = 3.88, 𝐿∗ = 6.4, 𝛥𝜌∗ = −0.01, 𝐺𝑎 = −955.4, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78, and 𝑡∗ = 1.
of the instabilities is lower. Fig. 6c shows the iso-dense case, which is
very similar to the ones shown in Fig. 4a and b, with an efficient and
stable displacement, and a thinner stagnant layer of the displaced fluid
near the walls and corners due to the level of stress, lower than the
7
yield stress. When the density difference is inverted, Fig. 6d, the stable
regime continues and fluid removal is nearly perfect, with an almost
negligible layer remaining adjacent to the walls. For this case, the final
efficiency is equal to 99.56%, as it can be observed in Fig. 7.
Journal of Non-Newtonian Fluid Mechanics 305 (2022) 104834P.R. Varges et al.
Fig. 13. (a) Velocity profile at the entrance of the enlarged region (𝑥∗ = 2𝑙∗ = 1.5) at
𝑡∗ = 1 for: 𝑅𝑒 = 3.88, 𝐿∗ = 6.4, 𝛥𝜌∗ = −0.01, 𝐺𝑎 = −955.4, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78.
From now on we present the effects of the other parameters to
better evaluate the role of the destabilizing mechanisms on the flow
displacement process.
4.2. Effect of Reynolds numbers
Fig. 8 shows the volume fraction after one pumped volume of fluid
2 at 𝑡∗ = 1, for 𝐿∗ = 6.4, 𝛥𝜌∗ = −0.01, 𝐺𝑎 = −955.4, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7,
𝑛 = 0.78 and: (a) 𝑅𝑒 = 0.04, (b) 𝑅𝑒 = 3.88, and (c) 𝑅𝑒 = 93.1, while
the displacement efficiency is depicted in Fig. 9. In addition, it also
shows the efficiency for 𝑅𝑒 = 3.88 and 𝑡∗ = 2 and 3 (blue squares).
The large value of 𝐺𝑎 indicates that buoyancy forces dominate over
viscous ones. For the lowest 𝑅𝑒, inertia is negligible and instabilities
are only due to buoyancy. It can be seen that the vortices are larger
than the ones shown in Fig. 6b, where the dimensionless yield stress
is higher. Therefore, viscous effects tend to dampen the instabilities
generated by buoyancy forces. As 𝑅𝑒 increases, 𝑅𝑒 = 3.88, inertia begins
to slightly affect, viscous dampening decreases and a small increase in
the displacement efficiency is observed. For the highest value of 𝑅𝑒,
instabilities are also observed, but the displacement is different due
to increased mixing between fluid, similar to what happened with the
iso-dense case (Fig. 4c). Inertia seems to overcome buoyancy, viscous
and yield stress effects, and promotes mixing between fluids so that the
displacing fluid removes fluid 1 more efficiently.
Fig. 10 shows the evolution in time of the volume fraction field for
the unstable case shown in Fig. 8b: 𝐿∗ = 6.4, 𝛥𝜌∗ = −0.01, 𝐺𝑎 = −955.4,
8
𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78 and 𝑅𝑒 = 3.88. Focusing the analysis in the
enlarged region, it is noted that the instability of the interface begins
as fluid 2 enters the enlarged region, with a pattern similar to the one
seen by Lee and Kim [33]. This interface evolves with time and remains
unstable. At the final time 𝑡∗ = 1, it is noted that a periodic pattern
of the interface is formed throughout the enlarged region. Patches of
fluid 1 (gray color) remain attached to the walls and fluid 2 (black
color) flows around them in a wavy path, establishing an unfavorable
displacement scenario.
To evaluate the effect of pumped volume, we analyzed the same
case, but increasing the pumped volume, up to three times the total
volume of the tube (𝑡∗ = 3), as depicted in Fig. 11. For 𝑡∗ = 2 it
is noted that near the wall at the exit of the enlarged region, fluid
2 become more concentrated due to buoyancy forces, which induces
fluid 1 to flow downstream, slightly increasing fluids mixing. This effect
continues as time evolves. Regarding the displacement efficiency, there
is an asymptotic improvement as time increases, due to the increased
mixing: it varies from 51.6% for 𝑡∗ = 1 to 58.9% for 𝑡∗ = 2 and 62.8%
for 𝑡∗ = 3, as can be seen in Fig. 9.
The velocity and rate-of-strain intensity fields at 𝑡∗ = 1 for this
same case are shown in Fig. 12. It is clear that the regions occupied by
fluid 1 are regions of almost zero strain rate (unyielded regions) and
velocity, indicating that fluid 1 remains stagnant. Looking at regions
occupied by fluid 2, there is a finger-type formation at the enlarged
region entrance. The velocity pattern of this finger repeats periodically.
Fig. 13 shows the velocity profile at the entrance of the enlarged region
(at 𝑥∗ = 2𝑙∗ = 1.5). It can be observed a parabolic shape for the
velocity profile though fluid 2 finger, close to the inner tube, with
a maximum deviated to the upper part of this finger. Moreover, the
velocities magnitude at the enlarged zone are very small, showing that
the fluid (mostly fluid 1 at this region) is almost stagnant. Accordingly,
we observe higher strain rates at the entrance of the enlarged region
and at the boundaries of the finger region.
4.3. Effect of viscosity ratio
Fig. 14 shows the volume fraction field at 𝑡∗ = 1 for 𝑅𝑒 = 3.88,
𝐿∗ = 6.4, 𝛥𝜌∗ = −0.01, 𝜏∗𝑦 = 12.7, and three viscosity ratios: (a)
𝜂∗ = 0.0008, (b) 𝜂∗ = 0.008, (c) 𝜂∗ = 0.08. It is important to observe
that in these cases, since the viscosity ratio varies, 𝑅𝑒 and 𝐺𝑎, which
were defined using fluid 2 properties, can give different insights when
we define them using the viscosity of fluid 1. So, if we look at these
quantities using the characteristic viscosity of fluid 1 instead of 𝜇2,
𝑅𝑒∗ = 𝑅𝑒.𝜂∗ and 𝐺𝑎∗ = 𝐺𝑎.𝜂∗. Therefore, 𝑅𝑒∗ = 3.1 × 10−3, 0.031
and 0.31, for cases (a), (b) and (c), respectively, and inertia effects
are low. At the same way, 𝐺𝑎∗ is equal to −1.3,−1.6 and −76.4, while
𝐺𝑎 = −1630,−196 and −955.4 for cases (a), (b) and (c), respectively.
For the lowest two viscosity ratios, the displacements are stable with
a large finger flowing through the enlarged region. Again, we observe
unyielded regions close to the walls and to the corners, where the shear
Fig. 14. Volume fraction at 𝑡∗ = 1 for 𝑅𝑒 = 3.88, 𝐿∗ = 6.4, 𝛥𝜌∗ = −0.01, 𝜏∗𝑦 = 12.7, and: (a) 𝜂∗ = 0.0008, (b) 𝜂∗ = 0.008, (c) 𝜂∗ = 0.08. Fluid 2 is represented by black color and fluid
1 by gray color.
Journal of Non-Newtonian Fluid Mechanics 305 (2022) 104834P.R. Varges et al.
Fig. 15. Displacement efficiency through time for 𝑅𝑒 = 3.88, 𝐿∗ = 6.4, 𝛥𝜌∗ = −0.01, 𝜏∗𝑦 = 12.7, and: (a) 𝜂∗ = 0.0008, (b) 𝜂∗ = 0.008, (c) 𝜂∗ = 0.08.
Fig. 16. Volume fraction at 𝑡∗ = 1 for 𝑅𝑒 = 3.88, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78, 𝛥𝜌∗ = −0.01, 𝐺𝑎 = −955.4, and : (a) 𝐿∗ = 1.8, (b) 𝐿∗ = 6.4, (c) 𝐿∗ = 13.7. Fluid 2 is represented by
black color and fluid 1 by gray color.
Fig. 17. Displacement efficiency through time for 𝑅𝑒 = 3.88, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 =
0.78, 𝛥𝜌∗ = −0.01, 𝐺𝑎 = −955, and : (a) 𝐿∗ = 1.8, (b) 𝐿∗ = 6.4, (c) 𝐿∗ = 13.7.
9
stress does not overcome the yield stress. It is interesting to observe that
the shape of the finger is slightly different between cases (a) and (b).
Inertia effects seems to be negligible in both fluids in all three cases.However, 𝐺𝑎∗ reveals that viscous forces in fluid 1 seem to balance
buoyancy for the first two cases, while for the case of 𝜂∗ = 0.08,
buoyancy forces are much stronger than viscous forces for fluid 1,
destabilizing the flow and reducing the displacement efficiency. The
resulting displacement efficiency through time is displayed in Fig. 15,
which shows the decrease of the efficiency as the viscosity ratio in-
creases. When instabilities are not present, the displacement efficiency
is very close for both cases, despite the fact of a 10 times increase in
the viscosity ratio. However, the onset of instabilities, generated by
buoyancy forces, causes a steep decrease in the displacement efficiency.
4.4. Effect of geometry
To analyze the effect of geometry on flow instabilities and displace-
ment efficiency, we present the volume fraction results for different
values of the dimensionless length 𝐿∗. Fig. 16 shows the volume
fraction field at 𝑡∗ = 1 for 𝑅𝑒 = 3.88, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78,
𝛥𝜌∗ = −0.01, 𝐺𝑎 = −955.4, and: (a) 𝐿∗ = 1.8, (b) 𝐿∗ = 6.4, (c)
𝐿∗ = 13.7. All cases are unstable flows due to buoyancy effects, with
similar patterns and periodicity, with a mild loose in periodicity as 𝐿∗
gets larger, probably due to a slight increase in mixing between fluids.
Therefore, the enlargement length does not affect flow stability, but
since mixing increases, the displacement efficiency improves, as it can
be seen in Fig. 17.
Journal of Non-Newtonian Fluid Mechanics 305 (2022) 104834P.R. Varges et al.
Fig. 18. Volume fraction at 𝑡∗ = 1 for 𝑅𝑒 = 3.88, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78, 𝐷∗ = 1.5, 𝐿∗ = 6.4, 𝛥𝜌∗ = −0.01, 𝐺𝑎 = −955.4, and: (a) 𝜎∗ = 0, (b) 𝜎∗ = 0.07, and (c) 𝜎∗ = 0.21. Fluid
2 is represented by black color and fluid 1 by gray color.
4.5. Effect of interfacial tension
Fig. 18 shows the effect of interfacial tension on flow stability. The
results are shown after one pumped volume of fluid 2, 𝑡∗ = 1, for 𝑅𝑒 =
3.88, 𝜂∗ = 0.08, 𝜏∗𝑦 = 12.7, 𝑛 = 0.78, 𝐿∗ = 6.4, 𝛥𝜌∗ = −0.01, 𝐺𝑎 = −955.4,
and 𝜎∗ = 0, 0.07, and 0.21. It is seen that the interfacial tension tends
to suppress the undulations and stabilize the flow, as expected. This is
more evident as flow develops upwards, and is stronger as interfacial
tension increases.
5. Conclusions
We analyzed numerically the displacement flow of a yield stress
material by a Newtonian one, through an enlarged region of vertical
annular tubes. This situation is usually found in cementing operations
of oil wells, where the drilling fluid (yield stress fluid) is pushed by a
spacer fluid (Newtonian) and this by the cement paste (also yield stress
fluid). The problem was solved using the finite volume technique and
the volume of fluid method to track the interface between phases.
The scaling showed that the flow can be described by several dimen-
sionless geometrical, flow and rheological parameters. Flow patterns
through the enlarged region, and the displacement efficiency were
obtained to investigate the role of the aforementioned parameters.
The results were obtained for typical properties, flow rates and
geometries of the oil industry, and mainly for situations where the
displaced fluid is heavier than the displacing one, representing the dis-
placement of the drilling fluid by the spacer fluid. Within this range of
parameters, instabilities, which are mostly caused by buoyancy forces,
have the greatest impact on the displacement process efficiency. In all
cases analyzed, the displacing fluid enters the enlarged region with a
finger-type configuration, and develops in different ways, depending on
the role of the buoyancy, inertia and viscous forces.
When viscous effects dominate, the finger-type configuration contin-
ues, and portions of the displaced fluid remain close to the walls, due
to yield stress effects. Flow patterns show that inertia and buoyancy
forces destabilize the flow, while viscous and interfacial forces tend
to damp these instabilities. The instabilities lead to a decrease on the
displacement efficiency. On the other hand, the dimensionless length
and inertia increase mixing, which in the end can be beneficial for the
displacement efficiency, although one has to have in mind that any
remaining displaced fluid can lead to contamination issues that may
impair the operation. This is important to point out, since there can
be situations where the displacement efficiency is high, but portions
of the displaced fluid remain mixed within the displacing fluid inside
the domain. It is worth recalling that any remaining fluid needs to be
washed later by the cement paste. Looking further for this stage of the
cementing operation, the wash out of the remaining fluid by the cement
paste will be easier when the remaining portions of the displaced
(drilling, yield stress) fluid are among the displacing one (spacer fluid)
than when the yield stress fluid remains close to the walls, because
stress levels are higher away from the walls. Therefore, comparing two
cases presenting highly efficient displacement, one obtained for higher
10
Reynolds number and another for lower Reynolds number, it would
probably be more interesting to work with highly inertia flows than
highly viscous ones.
The results obtained show interesting and important features of
the displacement of a yield stress material by a Newtonian one. The
development of unstable interfaces may lead to undesirable situations
that can result in an unpredictable and/or ineffective displacement
processes. The results have direct application in the oil industry, but
the analysis presented here is applicable for more general situations.
Declaration of competing interest
The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.
Acknowledgments
The authors are indebted to Petrobras, Brazil, CNPq, Brazil, CAPES,
Brazil, FAPERJ, Brazil, FINEP, Brazil, and MCT, Brazil for the financial
support to the Rheology Group at PUC-Rio.
Appendix A. Supplementary data
Supplementary material related to this article can be found online
at https://doi.org/10.1016/j.jnnfm.2022.104834.
References
[1] M.F. Naccache, H.A.M. Pinto, A. Abdu, Flow displacement in eroded regions
inside annular ducts, J. Braz. Soc. Mech. Sci. Eng. 40 (420) (2018) 1–14.
[2] A. Etrati, A. Roustaei, I.A. Frigaard, Strategies for mud-removal from washouts
during cementing of vertical surface casing, J. Pet. Sci. Eng. 195 (107454)
(2020).
[3] P.J.T. Espinoza, P.R. Varges, E.C. Rodrigues, M.F. Naccache, P.R. de Souza
Mendes, Displacement flow of yield stress materials in annular spaces of variable
cross section, J. Pet. Sci. Eng. 208 (109614) (2022).
[4] S. Taghavi, K. Alba, M. Moyers-Gonzalez, I. Frigaard, Incomplete fluid–fluid
displacement of yield stress fluids in near-horizontal pipes: Experiments and
theory, J. Non-Newton. Fluid Mech. 167–168 (2012) 59–74.
[5] K. Alba, S. Taghavi, J.R. de Bruyn, I. Frigaard, Incomplete fluid–fluid displace-
ment of yield-stress fluids. Part 2: Highly inclined pipes, J. Non-Newton. Fluid
Mech. 201 (2013) 80–93.
[6] T.G. Roberts, S.J. Cox, A.L. Lewis, S.A. Jones, Characterisation and optimisation
of foams for varicose vein sclerotherapy, Biorheology 57 (2–4) (2021) 77–85.
[7] P.R. de Souza Mendes, M.F. Naccache, P.R. Varges, F.H. Marchesini, Flow of vis-
coplastic liquids through axisymmetric expansions-contractions, J. Non-Newton.
Fluid Mech. 142 (1–3) (2007) 207–217.
[8] M.F. Naccache, R.S. Barbosa, Creeping flow of viscoplastic materials through a
planar expansion followed by a contraction, Mech. Res. Commun. 34 (5) (2007)
423–431.
[9] B. Nassar, P.R. de Souza Mendes, M.F. Naccache, Flow of elasto-viscoplastic
liquids through an axisymmetric expansion–contraction, J. Non-Newton. Fluid
Mech. 166 (2011) 386–394.
[10] P.R. Varges, B.S. Fonseca, P.R. de Souza Mendes, M.F. Naccache, C.R. Miranda,
Flow of yield stress materials through annularabrupt expansion–contractions,
Phys. Fluids 32 (083101) (2020) 1–12.
https://doi.org/10.1016/j.jnnfm.2022.104834
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb1
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb1
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb1
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb2
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb2
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb2
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb2
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb2
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb3
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb3
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb3
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb3
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb3
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb4
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb4
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb4
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb4
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb4
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb5
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb5
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb5
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb5
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb5
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb6
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb6
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb6
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb7
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb7
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb7
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb7
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb7
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb8
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb8
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb8
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb8
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb8
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb9
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb9
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb9
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb9
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb9
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb10
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb10
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb10
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb10
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb10
Journal of Non-Newtonian Fluid Mechanics 305 (2022) 104834P.R. Varges et al.
[11] A. Roustaei, A. Gosselin, I.A. Frigaard, Residual drilling mud during conditioning
of uneven boreholes in primary cementing. Part 1: Rheology and geometry effects
in non-inertial flows, J. Non-Newton. Fluid Mech. 220 (2015) 87–98.
[12] A. Roustaei, I.A. Frigaard, Residual drilling mud during conditioning of uneven
boreholes in primary cementing. Part 2: Steady laminar inertial flows, J.
Non-Newton. Fluid Mech. 226 (2015) 1–15.
[13] D.D. dos Santos, S.L. Frey, M.F. Naccache, P.R. de Souza Mendes, Flow of
elasto-viscoplastic liquids through a planar expansion–contraction, Rheol. Acta
53 (2014) 31–41.
[14] A. Roustaei, I. Frigaard, The occurrence of fouling layers in the flow of a yield
stress fluid along a wavy-walled channel, J. Non-Newton. Fluid Mech. 198 (2013)
109–124.
[15] W. Chin, X. Zhuang, Exact non-Newtonian flow analysis of yield stress fluids in
highly eccentric borehole annuli with pipe or casing translation and rotation,
in: SPE (Ed.), CPS/SPE International Oil & Gas Conference and Exhibition, no.
131234-PP in 1, 2010.
[16] R.C. Haut, R.J. Crook, Primary cementing: The mud displacement process, in:
SPE (Ed.), SPE, no. 8253 in 1, 1979.
[17] C.W. Sauer, Mud displacement during cementing: a state of art, J. Pet. Technol.
39 (1987) 1091–1101.
[18] C.F. Lockyear, A.P. Hibbert, Integrated primary cementing study defines key
factors for field success, J. Pet. Technol. 41 (1989) 1320–1325.
[19] T.R. Smith, K.M. Ravi, Investigation of drilling fluid properties to maximize
cement displacement efficiency, in: SPE (Ed.), SPE Annual Technical Conference
and Exhibition, no. 22775 in 1, 1991.
[20] T. Deawwanich, J. Liew, Q.D. Nguyen, M. Savery, N. Tonmukayakul, W. Chin,
Displacement of viscoplastic fluids in eccentric annuli: numerical simulation and
experimental vaidation, in: Chemeca 2008 Conference, Vol. 1, 2008, pp. 1–10.
[21] M. Savery, W. Chin, K.B. Yerubandi, Modeling cement placement using a new
3-D flow simulator, in: AADE (Ed.), AADE Fluids Conference and Exhibition, no.
AADE-08-DF-HO-08 in 1, 2008.
[22] P.E. Aranha, C.R. Miranda, J.V.M. Magalhaes, G. Campos, A.L. Martins, A.B.
Ramalho, M.F. Naccache, Dynamic aspects governing cement-plug placement in
deepwater wells, SPE Drill. Completion 1 (140144) (2011) 341–351.
[23] E.S.S. Dutra, M.F. Naccache, P.R. de Souza Mendes, A.L. Martins, C.R. Miranda,
Liquid displacement through tube-annular transition region inside oil wells, in:
ASME (Ed.), ASME International Mechanical Engineering Congress & Exposition,
no. IMECE2005-81279 in 1, 2005.
11
[24] A. Renteria, A. Maleki, I. Frigaard, B. Lund, A. Taghipour, J. Ytrehus, Effects
of irregularity on displacement flows in primary cementing of highly t deviated
wells, J. Pet. Sci. Eng. 172 (2019) 662–680.
[25] H.J. Skadsem, S. Kragseta, B. Lund, J.D. Ytrehus, A. Taghipour, Annular displace-
ment in a highly inclined irregular wellbore: Experimental and three-dimensional
numerical simulations, J. Pet. Sci. Eng. 172 (2019) 998–1013.
[26] S. Akbari, S.M. Taghavi, Fluid experiments on the dump bailing method in the
plug and abandonment of oil and gas wells, J. Pet. Sci. Eng. 205 (108920)
(2021).
[27] M.A. Tehrani, S.H. Bittleston, P.J.G. Long, Flow instabilities during annular
displacement of one non-Newtonian fluid by another, Exp. Fluids 14 (1993)
246–256.
[28] A. Renteria, P. Sarmadi, C. Thompson, I. Frigaard, Effects of wellbore irregularity
on primary cementing of horizontal wells, Part 1: Large scale effects, J. Pet. Sci.
Eng. 208 (109581) (2022).
[29] P. Sarmadi, A. Renteria, C. Thompson, I. Frigaard, Effects of wellbore irregularity
on primary cementing of horizontal wells, Part 2: Small scale effects, J. Pet. Sci.
Eng. 210 (110026) (2022).
[30] H. Jung, I.A. Frigaard, Evaluation of common cementing practices affecting
primary cementing quality, J. Pet. Sci. Eng. 208 (109622) (2022).
[31] D. Truzzolillo, L. Cipelletti, Hydrodynamic instabilities in miscible fluids, J. Phys.
Condens. Matter 30 (3) (2018) 033001.
[32] T. Atmakidis, E.Y. Kenig, A study on the Kelvin-Helmholtz instability using two
different computational fluid dynamics methods, J. Comput. Multiph. Flows 2
(1) (2010) 33–44.
[33] H.G. Leea, J. Kim, Two-dimensional Kelvin–Helmholtz instabilities of multi-
component fluids, Eur. J. Mech. B Fluids 49 (2015) 77–88.
[34] A.T. Dinh, H. Haraldsson, Z. Yang, B. Sehg, Advances in Fluid Mechanics III:
Simulation of Viscous Stabilization of Kelvin- Helmholtz Instability, WIT Press,
2000.
[35] C.W. Hirt, B.D. Nichols, Volume of fluid (vof) method for the dynamics of free
boundaries, J. Comput. Phys. 39 (1981) 204–225.
[36] Ansys Fluent Theory Guide, Release 2022 R1, Ansys Inc, 2022.
[37] Ansys Fluent User’s Guide, Release 2022 R1, Ansys Inc, 2022.
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb11
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb11
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb11
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb11
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb11
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb12
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb12
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb12
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb12
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb12
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb13
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb13http://refhub.elsevier.com/S0377-0257(22)00074-X/sb13
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb13
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb13
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb14
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb14
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb14
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb14
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb14
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb15
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb15
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb15
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb15
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb15
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb15
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb15
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb16
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb16
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb16
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb17
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb17
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb17
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb18
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb18
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb18
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb19
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb19
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb19
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb19
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb19
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb20
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb20
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb20
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb20
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb20
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb21
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb21
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb21
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb21
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb21
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb22
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb22
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb22
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb22
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb22
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb23
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb23
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb23
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb23
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb23
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb23
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb23
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb24
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb24
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb24
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb24
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb24
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb25
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb25
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb25
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb25
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb25
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb26
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb26
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb26
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb26
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb26
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb27
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb27
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb27
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb27
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb27
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb28
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb28
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb28
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb28
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb28
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb29
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb29
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb29
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb29
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb29
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb30
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb30
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb30
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb31
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb31
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb31
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb32
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb32
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb32
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb32
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb32
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb33
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb33
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb33
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb34
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb34
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb34
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb34
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb34
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb35
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb35
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb35
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb36
http://refhub.elsevier.com/S0377-0257(22)00074-X/sb37
	Flow instabilities in fluid displacement through enlarged regions in annular ducts
	Introduction
	Modeling
	Scaling
	Numerical solution
	Results and discussion
	Effect of density difference
	Effect of Reynolds numbers
	Effect of viscosity ratio
	Effect of geometry
	Effect of interfacial tension
	Conclusions
	Declaration of competing interest
	Acknowledgments
	Appendix A. Supplementary data
	References