J. Non-Newtonian Fluid Mech. 147 (2007) 92–108
Effect of contraction ratio upon viscoelastic flow
in contractions: The axisymmetric case
M
´
onica S.N. Oliveira
a
, Paulo J. Oliveira
b
, Fernando T. Pinho
c,d
, Manuel A. Alves
a,∗
a
Departamento de Engenharia Qu´ımica, CEFT, Faculdade de Engenharia da Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal
b
Departamento de Engenharia Electromecˆanica, Unidade de Materiais Tˆexteis e Papeleiros, Universidade da Beira Interior, 6201-001 Covilh˜a, Portugal
c
Centro de Estudos de Fen´omenos de Transporte, Faculdade de Engenharia da Universidade do Porto, 4200-465 Porto, Portugal
d
Universidade do Minho, Largo do Pa¸co, 4704-553 Braga, Portugal
Received 11 January 2007; received in revised form 10 July 2007; accepted 10 July 2007
Abstract
A comprehensive numerical study of the effects of the contraction ratio upon viscoelastic flow through axisymmetric contractions was carried out.
Six contraction ratios were examined (CR = 2, 4, 10, 20, 40 and 100) using the Oldroyd-B and Phan–Thien–Tanner (PTT) constitutive equations,
under creeping-flow conditions and for a wide range of Deborah numbers (De). The results enabled the construction of vortex pattern maps, with
CR and De as independent parameters, elucidating the role of these dimensionless groups in controlling vortex growth, vortex type (lip or corner
vortices), and pressure-drop characteristics. The extensional parameter of the PTT model was also varied (ε = 0–0.5) and it was found that for small
values of ε the Couette correction is a monotonic decreasing function of De, while for high ε values it is a monotonic increasing function.
© 2007 Elsevier B.V. All rights reserved.
Keywords: Axisymmetric contraction; Viscoelastic fluid; Creeping flow; Contraction ratio; PTT model; Oldroyd-B model
1. Introduction
The entry-flow is a long-standing problem dating back to
the late 1800s [1]. The early works were mainly experimental
and were primarily stimulated by the need to construct a capil-
lary rheometer capable of measuring accurately the viscosity of
Newtonian fluids [2]. Since then contraction flows have been the
subject matter of numerous experimental and numerical works.
The first numerical study of this problem, in which the complete
equations of motion were solved, was published by Vrentas et
al. [3] and concerned the creeping flow of a Newtonian fluid
through an axisymmetric contraction.
In spite of the simple geometry, these flows exhibit com-
plex patterns where shear and extensional regions co-exist. Near
the walls the flow is shear-dominated while along the center-
line it is purely extensional. These flows are amongst the most
studied extensionally dominated flows, since they assume partic-
ular importance in industrial applications involving viscoelastic
∗
Corresponding author. Fax: +351 225081449.
E-mail addresses: monica.oliveira@fe.up.pt (M.S.N. Oliveira), pjpo@ubi.pt
(P.J. Oliveira), fpinho@fe.up.pt, fpinho@dem.uminho.pt (F.T. Pinho),
mmalves@fe.up.pt (M.A. Alves).
non-Newtonian fluids. Examples worth mentioning are polymer
processing applications, such as injection molding, spinning and
film blowing [4]. Literature on contraction flows up until the
late 1980s have been the subject of detailed reviews by Boger
[1] and White et al. [5], in which the complex effects of flow
geometry and fluid rheology were addressed. For a comprehen-
sive account of more recent results the reader is referred to the
thorough introductions of Alves et al. [6] and Rodd et al. [7].
An overview of the evolution of numerical methods applied to
the flow of viscoelastic fluids through contractions can be found
in Walters and Webster [8], which includes typical constitutive
models used to represent real fluids; for a more complete review
on the theme refer to the book by Owens and Phillips [9]. Here
we limit ourselves to a description of some important results
with direct relevance to the problem at hand.
1.1. The sudden contraction configurations
There are three main types of configuration used in sud-
den contraction research: the axisymmetric, the planar and the
three-dimensional contraction of which the square–square con-
traction is a particular case. In the first, the fluid flows from a
capillary of large diameter through a sudden contraction into a
0377-0257/$ – see front matter © 2007 Elsevier B.V. All rights reserved.
doi:10.1016/j.jnnfm.2007.07.009
M.S.N. Oliveira et al. / J. Non-Newtonian Fluid Mech. 147 (2007) 92–108 93
smaller capillary; in the corresponding planar case, the capil-
laries are replaced by channels of large aspect ratio (quasi-2D)
and the contraction occurs only in one direction, while in the
square–square case this happens in two perpendicular directions.
The flow of viscoelastic fluids through sudden contractions,
either planar, axisymmetric or square–square, generates com-
plex flow patterns. In general, these flows comprise regions of
strong shearing close to the walls and non-homogeneous exten-
sion along the centerline upstream and downstream of the con-
traction [10]. Yet, it has been shown that the fluid behavior asso-
ciated with entry flows in different geometries can be quite con-
trasting as a result of geometric and rheological dissimilarities.
Perhaps the most widely studied configuration is the axisym-
metric contraction, mainly due to its implications in pipe/duct
flow [11]. Most of the early papers are devoted to the experimen-
tal study of this flow (e.g. [12–16]; more recent experimental
studies can be found in [10,17–20]). Up to the early 1990s,
numerical methods were usually unable to accurately predict the
steady contraction flow of viscoelastic fluids, and most efforts
were thus dedicated to improving the numerical techniques, as
documented in the books of Crochet et al. [21] and Owens and
Phillips [9].
Studies on the planar geometry were carried out soon after
the pioneering investigations on axisymmetric contractions. The
planar geometry assumed relevance as optical experimental
techniques evolved, such as flow visualization and especially
birefringence, and the number of published works on this config-
uration increased substantially [22–29]. The flow through planar
contractions has also received additional attention in numerical
investigations [30–39] since this configuration was chosen as
a test case for the assessment and improvement of numerical
methods in computational rheology [40].
Some of the findings observed for the axisymmetric geome-
tries are replicated in the planar case, such as the formation
of a recirculating region upstream of the contraction and the
existence of an extra pressure drop associated with the flow of
viscoelastic fluids. However, some major differences have been
identified as a result of the different strain and strain-rate his-
tories experienced by fluid elements in the two geometries [1].
Walters and Webster [41] found no significant vortex activity
for Boger fluids in the 4:1 planar case, in marked contrast to
observations in 4.4:1 circular contractions. However, for shear-
thinning fluids, vortex growth was observed in both planar and
axisymmetric geometries. Evans and Walters [25,26] studied
the flows of shear-thinning and constant-viscosity elastic fluids
and found strong vortex enhancement at all times for the shear-
thinning fluids, but once more no significant vortex activity was
observed for Boger fluids in planar contractions. Rothstein and
McKinley [10] found that in the axisymmetric case, the size of
the salient corner vortex formed was smaller than in the planar
case. They attributed this outcome to the different Hencky strains
experienced by the polymer molecules as the flow changes from
uniaxial to planar.
Over the years, little attention has been paid to the flow
through square–square contractions, most likely because they
are significantly more complex, and in numerical terms a “sim-
ple” 2D simulation is clearly inadequate. In many ways, the
square–square geometry can be thought of as an intermediate
between the axisymmetric and the planar contraction. Geomet-
rically, similarities between the planar and the square–square
contractions are easy to picture. On the other hand, in terms of the
actual flow and variation of strain-rates, similarities between cir-
cular and square–square contractions have been observed in the
experimental study of Walters and Rawlinson [42]. In their study,
the differences between the flow of Boger fluids in planar and
axisymmetric geometries were also seen to occur between pla-
nar and square–square contractions. More recently, a few studies
have addressed this geometry, both under an experimental [6]
and a numerical [43] perspective.
1.2. Flow patterns
One of the remarkable flow features of viscoelastic fluids
worth emphasizing is the vortex formation and vortex enhance-
ment mechanism upstream of the contraction plane. In general,
for axisymmetric contraction flows, strong vortex enhancement
is observed both for Boger fluids and for shear-thinning vis-
coelastic fluids. However, the vortex characteristics and the way
vortex enhancement evolves with the Deborah number can be
strikingly different, depending on the flow geometry and the fluid
rheology. In some cases the vortex forms near the salient corner
and increases in strength, growing steadily upstream and radi-
ally inward towards the re-entrant corner [1,18]. In other cases
the salient corner vortex grows in size with the Deborah num-
ber, while simultaneously a lip vortex forms near the re-entrant
corner; subsequently, the lip vortex grows radially outwards and
forces the corner vortex to decrease in both size and intensity
until it is completely overtaken by the lip vortex [10].
Understanding the mechanisms underlying the flow transi-
tions (and inherent vortex enhancement evolution) that take
place as the Deborah number is varied has been a main driv-
ing force behind the experimental work in this area [9]. This
evolution was shown to depend greatly on the fluid rheology.
Boger et al. [16] and Boger and Binnington [18] investigated
the behavior of two different Boger fluids having similar shear
properties and yet found quite different vortex dynamics. It was
then recognized that extensional properties had to be taken into
account, and Boger [1] suggested the primary parameter to be
the extensional viscosity. Recent experimental investigations by
Rothstein and McKinley [10,20] corroborate the important role
of the extensional viscosity on the dynamics of vortex growth
and associated enhanced pressure drop. In this context one of the
aims of the current paper is to examine under which conditions
one pathway of vortex evolution is preferred over the other.
1.3. Excess pressure drop
In addition to the kinematics, another essential flow charac-
teristic is the pressure drop resulting from a sudden change in
diameter. The total pressure drop is a result of the pressure drop
due to the fully developed viscous flow through the tubes plus
the excess pressure drop associated with contraction entrance
effects. Apart from being crucial for a proper assessment of the
pumping power required in ducts [44,45], it yields information
94 M.S.N. Oliveira et al. / J. Non-Newtonian Fluid Mech. 147 (2007) 92–108
about the global state of viscoelastic stresses in the flow [10].In
fact, many researchers have attempted to use the excess pressure
as a means to estimate the extensional viscosity of viscoelastic
fluids [19,46–51].
The extra pressure drop observed experimentally and pre-
dicted numerically for the flow of viscoelastic fluids through
sudden contractions can reach values much greater than those
exhibited by the corresponding Newtonian (or inelastic) fluids
with similar shear viscosity at the same flow rates [10,20,52].
This extra pressure drop is usually presented in dimensionless
form in terms of a Couette correction. For shear-thinning fluids,
the Couette correction is found numerically to attain a minimum
at low Deborah numbers and then to increase as the Deborah
number is further increased [37,38,53]. Furthermore, contrac-
tion ratio and geometry type have an effect on the magnitude and
on the onset of this enhanced pressure drop [10]. Experiments
with Boger fluids show a significant increase of the Couette
correction with the Deborah number, in marked contrast with
the numerical predictions using the UCM and the Oldroyd-
B models, for which a strong reduction in the extra pressure
drop is found, even leading to a pressure recovery due to elas-
tic effects [35,38]. Rothstein and McKinley [10,20] suggested
that the excess pressure drop observed in the experiments with
Boger fluids resulted from an extra dissipative contribution to
the elastic stress due to a stress-conformation hysteresis in the
non-homogeneous extensional flow near the contraction plane.
Phillips et al. [54] were able to predict numerically, at least qual-
itatively, a significant pressure drop enhancement using a closed
form of the adaptive length scale model of Ghosh et al. [55] in
accordance with experimental observations.
1.4. The contraction ratio
Ever since the fifth international workshop on numerical
methods in non-Newtonian flows [40], when the planar and
axisymmetric contraction geometries with a 4:1 contraction
ratio were put forth as benchmark problems in computational
rheology, special attention has been granted to geometries with
this specific contraction ratio. The reasoning behind the choice
of this particular value was based on the known results regarding
the flow of Newtonian fluids, where all interesting flow features
worthy of note take place for contraction ratios beneath 4:1
[52]. Even though it has been recognized a posteriori that for
non-Newtonian fluids this choice has been misguided, as many
interesting flow characteristics are only apparent for contraction
ratios beyond the specific ratio of 4:1 [52], most published
numerical works are still concerned with the standard case of
the 4:1 contraction ratio (e.g. [35,36,38,56–59]). Some authors
have used different values of the contraction ratio and there
are a few, mainly experimental, studies in which variations of
the contraction ratio are considered [1,10,15,17]. An exception
to this state of affairs in numerical investigations is the work
of Alves et al. [39], in which a detailed analysis of the effect
of the contraction ratio on the flow of a PTT fluid through a
planar sudden contraction was investigated. That work focused
on describing the flow patterns and quantifying the vortex
characteristics as a function of contraction ratio and Deborah
number. For high contraction ratios, it was found that flow
features in the vicinity of the corner, such as corner-vortex size,
and corresponding streamlines, scale with the upstream length
scale and with elasticity given by the Deborah number divided
by the contraction ratio, while features near the re-entrant corner
(vortex intensity and streamlines) scale with the downstream
length scale and with elasticity measured by the Deborah
number defined in terms of downstream quantities.
These experimental and numerical studies emphasize the
importance of investigating different contraction configurations
at varying contraction ratios. Thus, the present work extends
the previous study of Alves et al. [39], which was for the planar
geometry, by performing a comprehensive numerical analysis
on the effects of the contraction ratio on the flow patterns
and vortex dynamics in axisymmetric sudden contractions. A
comprehensive set of numerical simulations was performed,
ensuing new numerical results for the axisymmetric geometry
with the following contraction ratios: 2, 4, 10, 20, 40 and 100.
For the standard 4:1 ratio, the results given here for vortex
size and intensity have controlled accuracy and may therefore
be considered as benchmark data for the axisymmetric con-
figuration. A survey of literature for this case and comparison
of the existing data reveals as much scatter as that observed
for the planar case (shown in Fig. 1 of Alves et al. [38]). In
addition, the present work goes one step further by examining
other flow characteristics (such as excess pressure drop) and
by investigating various constitutive models corresponding
to Newtonian fluids, viscoelastic shear-thinning fluids (PTT
model) and constant-viscosity Boger fluids (Oldroyd-B model).
The paper is organized as follows. In Section 2, we present the
general governing equations and outline the numerical method
used to solve them. Section 3 describes the flow geometry and
the computational meshes used and in Section 4 we present
benchmark data for the standard 4:1 contraction geometry. In
Section 5, a systematic study of the effects of contraction ratio
and Deborah number on the flow patterns, vortex characteris-
tics and pressure drop is carried out, whereas Section 6 deals
with the effect of the rheological model. Finally, in Section 7,
we assess the suitability of the normal-stress ratio criterion of
Rothstein and McKinley [10] to estimate the onset of lip-vortex
activity, before ending the paper with the main conclusions.
2. Equations and solution method
The flow of a viscoelastic incompressible fluid can be
expressed in terms of its governing equations: conservation of
mass, conservation of momentum and a constitutive equation
for the extra stress tensor as a function of the flow kinematics.
The first two equations can be written as:
∇·u = 0 (1)
ρ
∂u
∂t
+∇·(uu)
=−∇p+η
s
∇·(∇u+∇u
T
)+∇ · (2)
where ρ is the density of the fluid, t the time, u the velocity
vector, p the pressure, η
s
the solvent viscosity and is the extra
stress tensor given by an appropriate constitutive equation.
M.S.N. Oliveira et al. / J. Non-Newtonian Fluid Mech. 147 (2007) 92–108 95
Throughout this paper, we will use three different constitutive
models: the Newtonian, the Oldroyd-B model and a simplified
form of the Phan–Thien–Tanner (PTT) constitutive equation
[60]. These three models can be written with the single following
equation:
λ
∂
∂t
+∇·(u)
+ f ()
= η
p
(∇u +∇u
T
) + λ( ·∇u +∇u
T
· ) (3)
where λ is the relaxation time of the fluid and η
p
is the polymer
zero-shear-rate viscosity (η
0
= η
p
+ η
s
, where η
0
is the total zero-
shear-rate viscosity).
The PTT model exhibits shear-thinning behavior and has
been used extensively in numerical studies involving contrac-
tion flows (see, e.g. [37–39] and references therein) and was
shown to be appropriate for modeling both polymeric solutions
and polymer melts [9], although often this may require the use
of a multimode version.
The simplified version used here does not include the
contribution from the lower-convected part of the full
Gordon–Schowalter derivative, and the linear form of the stress
function f() is selected,
f () = 1 +
λε
η
p
tr(), (4)
where tr( ) represents the trace operator and ε is an extension-
related parameter of the model, which eliminates singularities
in extensional viscosity. The limiting case of ε = 0 will also be
considered, as it represents an Oldroyd-B model, which is often
used to simulate the behavior of constant shear-viscosity Boger
fluids. Newtonian fluids are recovered by setting λ =0.
The governing equations, Eqs. (1)–(3), are solved using a
finite-volume method with a time-marching pressure-correction
algorithm with fully collocated variable arrangement, as detailed
in Oliveira et al. [61]. This methodology has been explained
thoroughly in previous works [35,38,39,62] and only a brief out-
line is given here. The governing equations are discretized on
a computational mesh, which in the present case is orthogonal
and non-uniform. The meshes used to map the computational
domain are created from structured blocks, which in turn are
organized in control-volumes or cells. The equations are then
discretized and transformed into systems of algebraic equations.
These relate the values of each dependent variable (p, u, ),
evaluated at the cell center, to the values in the neighboring
control-volume centers. The time-dependent terms in Eqs. (2)
and (3) are retained in the discretization so that a steady-state
solution is effectively approached by successive time advance-
ment steps. At each time step, the algebraic equations are solved
for the dependent variables, until the norm of the residuals of all
the equations drops below a convergence tolerance within the
range of 10
−4
to 10
−6
. In all cases full convergence was checked
by monitoring relevant variables, such as vortex size and inten-
sity, during the simulation procedure. Since in the present work
we are interested in steady-state calculations, the time derivative
is discretized with an implicit first-order Euler scheme. Central
differences are used to discretize the diffusive terms, while the
Fig. 1. Schematics of the flow configuration.
CUBISTA high-resolution scheme [62] is employed in the dis-
cretization of the advective term of the constitutive equation. The
advective term in the momentum equation in neglected since in
this work we are interested in truly creeping-flow conditions.
The CUBISTA scheme is a bounded version of the QUICK
scheme developed by Leonard [63], thus ensuring third-order
accuracy in uniform meshes, and was designed to yield a sym-
metric scheme and to respect Total Variation Diminishing (TVD)
restrictions. Alves et al. [39] used it extensively in a study of
PTT fluid flow through planar contractions, which is a similar
problem to the one we are dealing with in the present paper.
For the typical differential constitutive equations found in vis-
coelastic flows, which are hyperbolic in nature, the CUBISTA
scheme was shown to be more appropriate than other high-
resolution schemes from the point of view of accuracy and
iterative-convergence properties [62].
3. Problem definition and computational meshes
The geometry of the axisymmetric contraction is depicted in
Fig. 1, where R
1
and R
2
refer to the radii of the pipes upstream
and downstream of the contraction, respectively. The contraction
Table 1
Characteristics of the computational meshes
Mesh CR L
1
/R
2
L
2
/R
2
NC x
min
/R
2
=
r
min
/R
2
M1-CR2 2 2500 2500 7,732 0.020
M1-CR4 4 2500 2500 8,980 0.020
M3-CR4 4 2500 2500 35,920 0.010
M1-CR10 10 2500 2500 10,420 0.020
M1-CR20 20 2500 2500 11,956 0.020
M1-CR40 40 2500 2500 15,796 0.020
M1-CR100 100 5000 5000 23,920 0.020
M3-CR100 100 5000 5000 95,680 0.010
M1-CR4-SHORT
a
4 40 100 5,282 0.020
M2-CR4-SHORT
a
4 40 100 10,587 0.014
M3-CR4-SHORT
a
4 40 100 21,128 0.010
M4-CR4-SHORT
a
4 40 100 42,348 0.0071
M5-CR4-SHORT
a
4 40 100 84,512 0.0050
CR, contraction ratio; L
1
, length of upstream channel; L
2
, length of downstream
channel; NC, number of cells; x and r, cell spacing.
a
For CR = 4 very refined meshes were used in the benchmark simulations
with the Oldroyd-B fluid in order to assess the numerical uncertainty and the
true order of convergence. In these simulations smaller inlet and outlet tubes
were used.
96 M.S.N. Oliveira et al. / J. Non-Newtonian Fluid Mech. 147 (2007) 92–108
ratio between the two tubes, defined as CR = R
1
/R
2
, is an impor-
tant geometric parameter and in this study was varied between
2 and 100.
The total number of computational cells was varied depend-
ing on the contraction ratio used. Table 1 shows important
information on the computational domain and the meshes used,
such as the lengths of the upstream and downstream pipes (L
1
and L
2
), the total number of cells (NC) and the size of the
smallest cell normalized by the radius of the downstream pipe,
x
min
/R
2
= r
min
/R
2
, which is located at the re-entrant corner.
All meshes are orthogonal and non-uniform, with the size of
Fig. 2. Zoomed view of the axisymmetric computational meshes with contrac-
tion ratios (a) 4:1 (mesh M1-CR4), (b) 20:1 (mesh M1-CR20) and (c) 100:1
(mesh M1-CR100).
each cell relating to its neighbors by a geometric progression
within each direction, as illustrated in Fig. 2 for contraction
ratios of 4, 20 and 100 near the contraction plane. Cell cluster-
ing close to the re-entrant corner (shown in the zoomed views
of Fig. 2) and along the pipe walls was implemented since
the development of thin velocity and stress boundary layers is
expected to occur in these regions. The normalized minimum
cell spacing near the re-entrant corner is the same for all meshes
(x
min
/R
2
= r
min
/R
2
= 0.02) so that the local variations of the
resulting solution fields are not influenced by the mesh resolu-
tion. The exception is the benchmark 4:1 test case, for which we
have also used more refined meshes in order to assess the accu-
racy of the numerical results, and the CR = 100 case for which
a few tests were undertaken using mesh M3-CR100, in order to
further assess the numerical accuracy. The tube lengths L
1
and
L
2
were varied according to the contraction ratio to ensure fully
developed flow well upstream of the contraction plane and com-
plete flow redevelopment downstream of the contraction plane.
To ensure that this is achieved, we have chosen appropriate pipe
lengths based on our previous work for the planar case with a
4:1 contraction ratio [38].
The solvent viscosity ratio β = η
s
/η
0
will for most cases be
kept constant at the standard value of β = 1/9, while the exten-
sibility parameter of the PTT model will be usually taken as
ε = 0.25, which is typical of polymer melts or concentrated
solutions. PTT fluids exhibit shear-thinning behavior and an
extensional viscosity that is strain and strain-rate dependent.
The shear and extensional rheometrical properties for the PTT
model were presented in Fig. 23 of Ref. [38] and are not repeated
here for conciseness. In some cases the ε parameter will also
be varied, ranging from ε = 0.5 down to ε = 0, representing an
Oldroyd-B fluid.
The two other dimensionless groups characterizing the flow
are the Deborah and the Reynolds numbers, which are here
defined in terms of downstream characteristics,
De =
λU
2
R
2
(5)
Fig. 3. Estimated error in x
R
/R
2
as a function of mesh size for the Newtonian
fluid and for the Oldroyd-B fluid at De = 1 (meshes M1-CR4-SHORT to M5-
CR4-SHORT).
