Thickening of viscoelastic flow in a model porous medium

7.1 MB·

by  E. J. Hemingway

Thickening of viscoelastic flow in a model porous medium
E. J. Hemingway1, A. Clarke2, J. R. A. Pearson2 and S. M. Fielding1
1Department of Physics, Durham University, Science Laboratories, South Road, Durham, DH1 3LE, UK
2Schlumberger Gould Research, Madingley Road, Cambridge, UK, CB3 0EL, UK
(Dated: January 17, 2017)
Westudynumericallytwo-dimensionalcreepingviscoelasticflowpastabiperiodicsquarearrayof
cylinderswithintheOldroydB,fene-CRandfene-Pconstitutivemodelsofdilutepolymersolutions.
Ourresultscapturetheinitialmilddecreasethendramaticupturn(‘thickening’)seenexperimentally
inthedragcoefficientasafunctionofincreasingWeissenbergnumber. Bysystematicallyvaryingthe
porosityoftheflowgeometry,wedemonstratetwoqualitativelydifferentmechanismsunderpinning
this thickening effect: one that operates in the highly porous case of widely spaced obstacles, and
7 another for more densely packed obstacles, with a crossover between these two mechanisms at
1 intermediate porosities. We also briefly consider 2D creeping viscoelastic flow past a linear array
0 of cylinders confined to a channel. Our results here suggest that recent numerical reports [1] of
2 viscoelasticturbulenceinthis2Dgeometrymayneedtobecheckedforconvergencewithrespectto
n grid resolution, with a possible return to steady flow for finer grids.
a
J
I. INTRODUCTION ical value Wi . This upturn is often also accompanied
6 c
1 by the development of time-dependent flow fields, with
crossing streaklines and structure in the third spatial di-
] Flows of polymeric melts and solutions exhibit a rich mension [5, 7], into the page in the simplified sketch of
t phenomenology that has been widely studied [2]. Many
f Fig. 1. For low Weissenberg numbers a subtler effect is
o such materials exhibit viscoelastic properties which dra-
sometimes also seen, in which the pressure drop initially
s maticallyaffecttheirbehaviourinprocessesusedintheir
. decreases slightly relative to the Newtonian case [7] be-
t industrial application. Examples range from melt extru-
a fore the upturn just described for Wi>Wic.
m sion through coating process to flow in porous materi-
als. The latter example has particular relevance in the Some of these studies [11, 12] correlate the observed
-
d oil and gas industry where viscoelastic solutions are rou- upturn in the drag coefficient (described therein via an
n tinely employed for enhanced oil recovery, matrix stim- apparent viscosity) found for polymer solution flow in
o ulation and fracturing. For enhanced oil recovery, aque- outcrop rock samples with flows of the same solutions
c ouspolymersolutionshavelongbeenusedtocontrolthe observed in 2D microfluidic networks. These networks
[
Saffman-Taylorinstability[3,4]thatresultsfromthedis- comprise“pores”andrandomlysized“throats”onagrid

1 placement of viscous oil by lower viscosity brine within rotated 45◦ to the average flow direction. The experi-
v a porous reservoir rock. Since the earliest experiments ments were subsequently extended [8] to show how the
3 it has been known that for certain polymeric solutions, onset of additional drag depended on solution parame-
3
even for single-phase flow, an anomalously high pressure ters. In particular, the onset of thickening was found to
2
gradient is observed when compared with the flow of an be well characterized by a Weissenberg number derived
4
0 appropriateequivalentNewtoniansolution. Whereasthis usingacharacteristicapparentshearrateproportionalto
. phenomenon has been long known, no clear understand- flow rate together with a molecular relaxation time (see
1
ing of the root cause of the high observed drag has been Fig. 11 in that work). In these studies the transition to
0
7 forthcoming. In particular there has been no clear eluci- a time-dependent, turbulent state is marked in the mi-
1 dation of the relative contributions of shear, extensional crofluidic networks by the appearance of crossing streak-
: andelasticstresses. Theexcesspressuregradientdirectly lines, a 3D effect which isn’t possible in the current 2D
v
impacts the industrial process by limiting the rate of in- study. NMR studies [13] demonstrated time-dependent
i
X jectionofpolymericsolution,whereanexcessivepressure flows within a full 3D pore network (a rock) via an effec-
r gradientwillleadtofracturingoftherockinthevicinity tivediffusionconstantmeasurement. Theanalysisdevel-
a of the injector and to degradation of the polymer. oped in the following paper can also be applied to these
geometries at 45 degree orientations. Our results (not
Experimentally, several authors have studied the flow
discussed here) suggest that the cylinder size has less ef-
of a viscoelastic fluid past a biperiodic array of obstacles
fect on the flow character q (defined in Section VIA)
arranged in two spatial dimensions [5–10]. (The upper
than for 0◦ orientations. However we defer a full study
panel of Fig. 1 shows a sketch of the simplified cartoon
in these rotated geometries to future work.
of such a geometry that we shall study numerically in
thiswork.) Thedominantphysicaleffectconsistentlyre- Fromatheoreticalviewpoint,earlyattemptstounder-
ported in these experiments is a dramatic upturn in the stand flow in porous media adopted a coarse-grained ap-
adimensional pressure drop, known as the drag coeffi- proach,discardingmicroscopicdetailsinfavourofmacro-
cient, relative to that for a Newtonian fluid of matched scopic properties. For Newtonian flow, Darcy [14] pro-
viscosity, for Weissenberg numbers Wi exceeding a crit- posed a relation between the pressure drop per unit
2
length of material ∆P/L and the mean velocity scale V nificant efforts have also been devoted to understanding
viscoelastic flow past a single cylinder or linear array of
∆P ηV cylinders confined to a channel: experimentally [19, 25–
= , (1)
L K 28], by linear stability analysis [24, 29], and by direct
numerical simulation [1, 19, 27, 30–34]. The lower panel
where η is the fluid’s viscosity. The permeability K is
of Fig. 1 shows a sketch of the simplified cartoon of such
a constant that should depend only on the properties of
a geometry that we shall study numerically in this work.
themedium,butisunknownapriori. Arelationbetween
BystudyingtheflowofBogerfluidspastasinglecylin-
the permeability and the medium’s porosity (cid:15) (the ratio
der in a channel, McKinley et al. [25] observed a transi-
of free volume to total volume) was later proposed in
tionfromsteady2Dtosteady3Dflowinthedownstream
the Blake-Kozney-Carman equation [15], which proved
wake. At higher flow rates, they found another insta-
successful in describing a range of simple flows. However
bility where time-dependent velocity oscillations form in
thevalidityofthismacroscopicapproachremainslargely
the wake region. Liu [19] considered flow past a single
limited to Newtonian flows.
cylinder, and widely and closely spaced linear arrays of
In theoretically understanding viscoelastic creeping
cylinders. For the single cylinder, they observed a mild
flows in porous media, much of the progress has been
downturn in the drag at moderate Wi, followed by an
made computationally [16–22]. Early simulations of
upturn at larger Wi which was accompanied by a tran-
two-dimensional (2D) viscoelastic flow past a bi-periodic
sition from steady 2D to steady 3D flow. In the lin-
square [16, 17] or rectangular [18] array of cylinders ob-
ear arrays, the transition was from steady 2D to time-
served a slightly reduced pressure drop at low Weis-
dependent 3D structure in both cases. Moss and Roth-
senberg numbers compared to that of a Newtonian fluid
stein[26]studiedflowofwormlike-micelles(WLMs)past
of matched viscosity, as seen experimentally. However
asinglecylinderforseveralratiosofcylinderdiameterto
thedramaticupturninthedragcoefficientathighWi>
channelwidth. Theyobserveasignificantdecreaseinthe
Wi , which is the dominant physical effect seen exper-
c
normalised pressure drop as a function of Wi which was
imentally, was not captured in these early numerical
attributed to the shear-thinning properties of the fluids.
works, presumably due to the restricted computational
Of the two fluids tested, only one exhibited an instabil-
processing power available at the time. It was however
ity, which was attributed to breakdown of the WLMs in
later captured in simulations of the Oldroyd B and fene-
the extensional flow in the cylinder wake. Using flow-
CR models, also in 2D biperiodic arrays [19].
inducedbirefringencemeasurements,theauthorsshowed
Alcocer et al. [20, 21] investigated 2D flow of the fene-
that shear flows at the channel walls were not necessary
CR model past a biperiodic array of cylinders, with a
to produce the wake instability. Recent similar experi-
particular focus on the dependence of the effective per-
ments [28] found that upstream vortices formed at much
meability on the cell aspect ratio, for a fixed area frac-
larger Wi∼103, and unsteady flow downstream at even
tion of cylinders. They demonstrated a non-monotonic
larger Wi∼104.
dependence of permeability on aspect ratio. For a fixed
aspectratioandareafraction,theyreportedinaninitial Simulationsof2Dcreepingviscoelasticflowpastalin-
increase in permeability with increasing Wi, equivalent ear array of cylinders confined in a channel [30, 31] cap-
to the initial decrease in drag in other studies. tured an initial mild decrease then subsequent upturn
Gillissen [22] simulated 2D flow of the fene-P model in the drag coefficient as a function of increasing Weis-
past a biperiodic hexagonal array of cylinders. This senberg number. All the states observed were however
study convincingly captured both the initial downturn steady as a function of time. In 2D simulations of flow
and then significant upturn in the drag coefficient seen pastasinglecylinderconfinedinachannel, temporalos-
experimentally as a function of Weissenberg number. cillations in the size of a recirculating region that forms
Analysing the flow field as a function of space in terms downstream of the cylinder are seen [32, 33], with an as-
ofregionsofpureshear(whichisthesameasextension), sociated slight increase in the drag coefficient compared
simpleshearandpurerotation,theydemonstratedapre- with the time-independent state.
dominanceofshearregionsatlowWi,withaprogressive AttemptstounderstandtheonsetathighWiofthe3D
increase in elongational regions with increasing Wi. At time-dependent states seen experimentally for creeping
high Wi the polymer conformation tensor was found to viscoelastic flow past an array of cylinders in a channel
befullyextended,showingtheimportanceoffinitechain have been made by performing a linear stability anal-
extensibility in this regime. ysis for the dynamics of small amplitude 3D perturba-
De et al. studied 3D flow of a FENE-P fluid past an tions to an initially 2D flow state. By such an analysis,
array of cylinders (both with and without walls) [23]. Smith et al. [24] reproduced some of the experimentally
They found an elastic instability whereby recirculating observed characteristics of the instability, in particular
regions in the cylinder wake break symmetry and form a the wavevector of the most unstable mode and the crit-
3D structure, which occurs at a Deborah number consis- ical Wi at which the instability first arises. Sahin and
tent with Ref. [24]. All their simulation runs attained a Wilson [29] demonstrated that the wavelength of the in-
time-independent steady-state stability scales with cylinder spacing for closely spaced
Besides the biperiodic geometries just discussed, sig- cylinders,andwiththesizeofthewakebehindthecylin-
3
der for wider cylinder spacing. Both studies warn that outline our numerical methods, and provide benchmarks
this 3D instability potentially restricts the range of Wi tovalidatethemagainstknownresultsintheNewtonian
over which purely 2D simulations might remain valid. limit. Wethenpresentourresultsforviscoelasticflow: in
V´azquez-Quesada and Ellero performed 2D smoothed Sec.VIforabiperiodicarrayofcylindersandinSec.VII
particle hydrodynamics (SPH) simulations of viscoelas- for a linear array of cylinders in a channel.
tic flow past an array of cylinders in a channel [1, 34],
capturingtheinitialdownturnandsubsequentupturnin
thedragcoefficientasafunctionofWeissenbergnumber, II. FLOW GEOMETRIES
as seen experimentally. At higher Wi the results of this
studydepartedfromother2Dnumericalworks[19]inre- Weshallstudytwodifferentflowgeometries. Thefirst
portingatransitiontoatime-dependentstate,whichthe comprisesabiperiodicarrayofcylinders, sketchedinthe
authors interpreted as viscoelastic turbulence. We shall upper panel of Fig. 1. The second comprises a linear
returntowardstheendofthispapertoappraisethisfind- array of cylinders in a channel bounded by solid walls,
ingusingourownnumericalmethod. Ourresultssuggest sketched in the lower panel of the same figure. In each
thatanyinstabilitytoatime-dependentstatemaydisap- case the flow cell has length L horizontally and height
x
pear with increasing mesh refinement in 2D simulations. L vertically, and the cylindrical obstacle has radius R.
y
Ribeiroetal.[27]studiedfully3Dviscoelasticflowpast Attheboundariesofthecellsrepresentedbydashedlines
acylinderconfinedtoachannel,bothexperimentallyand all flow variables are ascribed periodic boundary condi-
numerically, for both a shear thinning and a Boger fluid. tions. At the cell boundaries represented by solid lines,
For the shear-thinning fluid they reported an elastic in- and at the cylinder surface, conditions of no-slip and no-
stability setting in upstream of the cylinder at critical permeation apply. In each case we assume the flow to
value of Wi that depends on the cylinder height in the be translationally invariant into the page in Fig. 1, and
vorticity direction. For Weissenberg numbers just be- simulatetheflowinthetwodimensionsofthepageonly.
yond onset of the instability, the system’s state is asym- Theflow,whichweassumetobefromlefttoright,can
metricandtime-independent. Asubsequenttransitionto be imposed in two different ways. In the first, a given
a time-dependent state was reported at larger Wi still. throughput Q per unit time is prescribed (such that the
characteristic velocity scale is then V =Q/L ), with the
Among the simulation studies of viscoelastic flow just y
pressure drop ∆P measured in response. Alternatively,
surveyed [16–22], each considered one (or in some cases
we can impose the pressure drop ∆P and measure the
two) fixed ratio(s) of obstacle size to obstacle spacing,
resulting throughput Q.
i.e., a fixed medium porosity. A key contribution of this
In either case, a key experimental observable is the
work is systematically to vary the medium porosity over
drag coefficient
a broad range, from the limit of widely spaced obsta-
cles to ones that nearly touch. In doing so, we shall ∆P
C = , (2)
demonstrate two qualitatively different mechanisms un- D ηV
derpinning the dramatic upturn of the drag coefficient
which measures the pressure drop normalised by the
with Weissenberg number seen experimentally: one at
characteristic velocity scale and the solvent viscosity. In
low obstacle area fraction and another at high area frac-
thelimitofNewtonianflowthisquantitydependsonlyon
tion, with a crossover between these mechanisms at in-
theflowgeometry. Innon-Newtonianflowitalsodepends
termediate area fraction.
via the Weissenberg number on the nonlinear constitu-
In neither case, however, do we find the thickening to
tive behaviour of the fluid in question. In reporting our
beassociatedwiththeonsetofatime-dependentflow,as
numerical results below we shall typically show the drag
often reported experimentally. The same is true in the
coefficient at any given Weissenberg number, C (Wi),
2Dbiperiodicsimulationsintheearlierliterature[16–22]. D
normalised by the corresponding value in the Newtonian
The most likely explanation for this discrepancy is that
limitofzeroWeissenbergnumber,C (Wi→0),defining
thetime-dependentstateseenexperimentallyis3Dinna- D
ture, and cannot be captured in purely 2D simulations. C (Wi)
χ= D . (3)
As noted above, in 2D simulations of flow past a linear C (Wi→0)
D
array of cylinders bounded by walls, V´azquez-Quesada
and Ellero [1, 34] did recently report a transition to a
time-dependent and apparently turbulent flow at high
III. CONSTITUTIVE MODELS
Weissenberg number. An additional contribution of our
work is to suggest that this time-dependence may disap-
We write the total stress T(r,t) in a fluid element at
pear either with (i) sufficient numerical mesh refinement
position r and time t as the sum of a viscoelastic con-
or (ii) by strictly enforcing fluid incompressibility.
tribution Σ(r,t) from the polymer chains, a Newtonian
The paper is structured as follows. In Secs. II and III
solvent contribution of viscosity η, and an isotropic con-
we introduce the flow geometries and constitutive mod-
tribution with a pressure p(r,t):
els to be studied. The governing parameters and dimen-
sionless groups are summarised in Sec. IV. In Sec. V we T =Σ+2ηD−pI. (4)
4
with a constant modulus G. In addition to the spring
force, each bead also experiences viscous drag against
thesolvent[35]andstochasticthermalfluctuations. The
conformation tensor is then taken to obey
∇ 1 (cid:96)2
W =− [f(W)W −g(W)I]+ ∇2W, (9)
τ τ
with a characteristic relaxation time τ, where
∇
W ≡(∂ +v·∇)W −KT ·W −W ·K (10)
t
is the upper convected derivative [35]. We have included
a modification to the original equations by introducing a
diffusiveterm,where(cid:96)isasmalllengthscalebelowwhich
gradients in W are attenuated. Similar modifications
have been made to the Johnson-Segalman model in the
contextofshearbanding[36,37]. Specificallyinthecon-
textofporousmedia,Gillissenincludeddiffusivetermsin
hisstudyofaFENEfluidpastabiperiodicarrayofcylin-
ders [22]. Thomases et al. also recently showed that a
small diffusive contribution can support a finite polymer
stress in a qualitatively similar fashion to FENE models
[38, 39]. Strictly, gradient terms in W require a bound-
FIG. 1. Flow geometries to be studied. Upper: Biperiodic aryconditionatthewallsandatthecylindersurface. For
array of cylinders. Lower: Linear array of cylinders in a simplicity we choose zero gradient at the walls ∂yW =0
channel. Solid lines represent closed walls, and dashed lines (when present). We fix (cid:96) = 0.01 in all that follows and
represent periodic boundaries. The flow is assumed transla- we have checked that our results are qualitatively (and
tionally invariant into the page. inalmostallcasesquantitatively1)unchangedbyfurther
decreasing (cid:96).
The functions f(W) and g(W) in Eqns. 8 and 9 are
ThesymmetricstrainratetensorD = 1(K+KT)where
2 different in the three different models. The Oldroyd B
K =∂ v and v(r,t) is the fluid velocity field.
αβ α β model has f(W) = g(W) = 1, which corresponds to
Throughoutweconsiderthecreepingflowlimitofzero
assuming that the spring of each dumbbell is Hookean.
Reynolds number. Here the condition of force balance
Forasustainedimposedextensionalstrainrate(cid:15)˙>1/2τ
requires the stress field T(r,t) to be divergence free:
this model displays an extensional catastrophe in which
thedumbbellsstretchoutindefinitelyandtheextensional
∇· T =0, (5)
stress grows indefinitely.
such that Thephenomenologicalfene-CRandfene-Pmodelsreg-
ularise this catastrophe by insisting that the extensional
∇.Σ+η∇2v−∇p=0. (6)
stress of the polymer chains (dumbbells) must remain fi-
nite at all deformation rates. They do so by replacing
The pressure field p(r,t) is determined by enforcing flow
the Hookean spring law by a non-linear law with finite-
incompressibility:
extensibility [40]. The fene-CR model has
∇·v =0. (7)
1
f(W)=g(W)=α(W)≡ , (11)
ThedynamicsofthepolymericstressΣisspecifiedby 1− Tr(W)
a viscoelastic constitutive model. In this work we con- Λ2
siderthreedifferentphenomenologicalconstitutiveequa- while the fene-P model has
tions: the Oldroyd B, fene-CR and fene-P models [35].
f(W)=α(W) and g(W)=1. (12)
These each describe a dilute polymer solution by rep-
resenting each polymer chain as a simplified dumbbell
The parameter Λ in Eqn. 11 characterises the max-
comprising two beads connected by a spring. The con-
imum extent to which any dumbbell can be stretched.
formation tensor W = (cid:104)RR(cid:105) is defined as the ensemble
average (cid:104)(cid:105) of the outer dyad of the dumbbell end-to-end
vector R, which is taken to have unit length in equilib-
rium. The conformation of the polymer chains deter- 1 Theonlyexceptionisinthecloselyspacedcylinderscaseshown
mines the viscoelastic stress according to inFig.15(right)wherethedragupturnisslightlylesssteepfor
(cid:96)=0.005,0.0025. Thepointofupturninthedragisunchanged
Σ=G[f(W)W −g(W)I], (8) inallcases,andallstatesremaintime-independent.
5
functionofγ˙. Thefene-Pmodelshearthinstoanextent
determined by the value of δ.
2
10 Inanidealviscometricplanarextensionalflow,theim-
posed velocity gradient tensor
/G
σ
(cid:18) (cid:19)
1 0
∇v| =(cid:15)˙ . (14)
ext 0 −1
0
10
In steady state, the extensional stress Σ −Σ (which
xx yy
we also denote σ ) as a function of (cid:15)˙ defines the fluid’s
E
extensional constitutive curve, as shown in Fig. 2 for the
-2 three models. As can be seen, the extensional constitu-
10
-2 -1 0 . 1 2 3
10 10 10 10 10 10 tive curve of the Oldroyd B model is undefined for flow
γ τ
rates (cid:15)˙ > 1/2τ, consistent with our discussion above of
thechainstretchcatastrophe. Incontrast,thefenemod-
104 elshavewell-definedconstitutivecurvesatallstrainrates
(cid:15)˙. Indeed for matched δ the two fene models have the
same extensional constitutive curves, with a finite limit-
ing extensional viscosity 2Gτ/δ as (cid:15)˙→∞.
2
10
/G
σ
IV. PARAMETERS AND DIMENSIONLESS
e
0 GROUPS
10
The eight parameters characterising the fluid, geome-
try, andimposedflowjustdescribedare: thesolventvis-
-2
10
-2 -1 . 0 1 cosity η, the polymer modulus G, the polymer viscoelas-
10 10 10 10
ε τ ticrelaxationtimescaleτ,thepolymerfiniteextensibility
parameter δ, the cell length L , the cell height L , the
x y
cylinder radius R and the flow’s throughput rate Q.
FIG. 2. Upper: Stationary constitutive curves for homoge-
neous simple shear flow in the Oldroyd B model (solid black We are free to choose units of mass, length and time,
curve), the fene-CR for δ = 0.001 (magenta dotted curve) leaving five dimensionless groups as follows: the ratio of
and 0.01 (magenta dashed curve) and the fene-P model for cylinder radius to gap height R˜ =R/L , the ratio of cell
y
δ = 0.001 (green dotted curve) and δ = 0.01 (green dashed length to cell height L˜ = L /L , the ratio of solvent
x x y
curve). The fene-CR curves are indistinguishable from the
to (zero shear) polymer viscosity β = η/Gτ, the finite
Oldroyd B curve.
extensibilityparameterδ,andaWeissenbergnumberWi
Lower: Counterpart stationary constitutive curves for ho-
characterisingthestrengthofthevelocitygradientscom-
mogeneous planar extensional flow, for the same parameter
pared to the inverse of the fluid’s stress relaxation time
values as in the upper panel, and with the same line key. In
τ. (We return below to define Wi precisely.) We drop
thiscase,formatchedδ,thefene-CRandfene-Pcurvesarein-
distinguishablefromeachother(andthetwomagentacurves tildes hereafter with the understanding that (for exam-
are accordingly hidden by the green ones). ple)valuesforRandLx arealwaysquotedinunitsofLy.
To allow benchmarking of our results against the earlier
literature we fix the viscosity ratio β = 0.59 through-
We choose to express this in terms of δ = Λ−2, noting out. Remaining to be explored numerically are then the
thatthelimitδ →0correspondstoOldroydBdynamics dimensionless cylinder radius R, the dimensionless cell
with infinite extensibility. length Lx, the finite extensibility parameter δ and the
Underconditionsofidealviscometricsimpleshearflow, Weissenberg number Wi.
theimposedvelocitygradienttensorinCartesian(x−y) The Weissenberg number is a dimensionless quantity
coordinates is characterising the scale of velocity gradients in the flow
inunitsoftheinverserelaxationtime. Onpurelydimen-
(cid:18) (cid:19)
0 0 sional grounds, one simple possible definition is Vτ/L .
∇v| =γ˙ . (13) y
shear 1 0 However we have found it more useful to define two dif-
ferent Weissenberg numbers based on a characterisation
In steady state under this applied flow, the shear stress of the velocity gradients that actually develop inside the
Σ (which we also denote σ) as a function of γ˙ defines flowgeometry. WeshallreturntodiscusstheseinSec.VI
xy
thefluid’sshearconstitutivecurve,asshowninFig.2for below.
the three models considered here. The shear viscosity In our studies of the biperiodic geometry we consider
of the Oldroyd B and fene-CR models is constant as a a square cell with L = 1.0. For the channel geometry
x
6
3
wetakethecylinderspacingL asanadditionalvariable 10
x
to be explored numerically.
V. NUMERICAL METHODS
2
C 10
Wenumericallysolvethemodelequations6,7,8and9 D
usingatimesteppingapproach. Eachmaintimestepcom-
prises two separate substeps, as follows. In the first sub-
step the viscoelastic stress Σ(r,t) is updated according
to the viscoelastic constitutive equations 8 and 9 with a
fixed velocity field v(r,t). This update is performed in 101
thewholeplaneofFig.1,eveninsidethecylinder,where
0 0.1 0.2 0.3 0.4 0.5
the velocity is zero to within the accuracy of our numer- R
ical approach. For this substep we adopt a method that 3
we have used in our own previous work [41], to which we
refer the reader for details.
With that newly updated viscoelastic stress field, the 2
velocity field is then updated in the second substep to
satisfy Eqns. 6 and 7. This substep needs more careful
discussion, because the presence of the cylindrical posts
1
renderstheproceduremorecomplicatedthaninRef.[41]
(which considered a rectangular flow cell with no obsta-
cles). In particular, the boundary conditions of no-slip 0
and no-permeation must be satisfied for the fluid veloc-
ity round the edge of each obstacle. To tackle this we FIG. 3. Upper: Drag coefficient as a function of cylinder
use an immersed boundary method (IBM)[42–44], which radiusforNewtonianflowinthebiperiodicgeometry. Dotted
couples a solution of Eqns. 6 and 7 on a regular Eulerian line: analytical prediction in the limit R → 0 [45]. Dashed
line: analytical prediction in limit R → L of near-touching
grid with an off-grid Lagrangian description of the cylin- y
cylinders [45]. Symbols: numerical results using immersed
der surface. (We use the term Lagrangian to be consis-
boundary method (circles), propagator method (squares),
tent with the IBM literature, although in our particular
phase field method with l = 0.02 (diamonds) and l = 0.01
problem the cylinder doesn’t convect with the flow.)
(triangles). (See [46] for definition of the parameter l.)
The cylinder surface is characterised by a one- Lower: Streamlines (red lines) and map of velocity magni-
dimensional curvilinear Lagrangian coordinate ξ around tude |v| calculated using IBM (left) and propagator (right)
it. (Although we use the word cylinder, in our 2D study methods for a cylinder radius R=0.25 and overall through-
itisofcourserepresentedbyacircularcrosssection,with put rate Q=1.0.
a 1D edge.) The location of a point at the Lagrangian
coordinate ξ is given in the Eulerian frame as X(ξ,t).
A Lagrangian force density F(ξ,t) is then incorporated,
then gives a force contribution in the Eulerian frame of
calculatedbyprescribingthedesiredlocationofthecylin-
derX (ξ),andimposingaHookeanrestoringforcewhen
0 (cid:90)
the cylinder deviates from this:
f(r,t)= F(ξ,t)δ(r−X(ξ,t))dξ. (18)
Ω
F(ξ,t)=−κ(X(ξ,t)−X (ξ)), (15)
0
where κ is a large spring constant. Given that the de- This is incorporated as an additional source term to the
sired location X is independent of time, differentiating lefthandsideofthegeneralisedStokes’equation6,which
0
Eqn 15 with respect to time gives is then solved on a rectangular Eulerian grid in the full
plane of Fig. 1 using the methods of Ref. [41]. With-
∂F(ξ,t)
out walls, this can be done either imposing the overall
=−κU(ξ,t), (16)
∂t throughout(andmeasuringthepressuredrop)orimpos-
in which the Lagrangian velocity U at the cylinder sur- ing the pressure drop (and measuring the throughput),
faceiscalculatedfromthevelocityv(r,t)intheEulerian and we have checked that these are equivalent in all our
frame as biperiodic results presented below. With walls included,
we always impose the latter quantity.
(cid:90)
U(ξ,t)= v(r,t)δ(r−X(ξ,t))dx. (17) The transfer of information between the Lagrangian
Ω and Eulerian grids in Eqns. 17 and 18 is achieved in nu-
Eqn. 16 is evolved at each timestep using an explicit merical practice by approximating the Dirac delta func-
Euler algorithm. The Lagrangian force density F(ξ,t) tion δ(r) by a smoothed Peskin delta function δ (r) =
P
7
δx(x)δy(y), in which [47]
P P numerical (L = 4)
x
(cid:115) 102 analytical, Faxen (1946)
|r| |r| (cid:18)|r|(cid:19)2
8hδ (r)=3−2 + +1+4 −4 |r|≤h
P h h h
(cid:115) C
|r| |r| (cid:18)|r|(cid:19)2 D
=5−2 − −7+12 −4 h≤|r|≤2h
h h h
=0 otherwise.
1
10
Here h = ∆x = L /N = ∆y = L /N is the Eulerian
x x y y
-2 -1
gridspacing,givenarectangulargridof(Nx,Ny)points. 10 R / L 10
The Lagrangian boundary of circumference 2πR is dis-
cretizedintoM nodesofequalseparation∆s=2πR/M.
Theoptimalvalueoftheratioα=∆s/hisunknownapri-
ori. Too small a value risks oversampling the boundary
forces, while too large a value risks fluid leakage across
the boundary. We set α=2, and have checked that our
results are robust to reasonable variations around this
value. We have also ensured that all our results pre-
sented below are converged on the limit of grid spacing
h→0 and of timestep ∆t→0.
In the numerical solution just described, the methods
employedtoupdatetheviscoelasticstressinthefirstsub-
stepofeachtimestephavebeentestedandbenchmarked
by ourselves in several previous publications [41, 50].
Therefore we focus here only on benchmarks to validate
the second substep, in which the velocity field is calcu-
lated by solving Eqns. 6 and 7. The presence of vis-
coelasticity in this substep is trivial: it appears only as FIG. 4. Upper: Comparison of the drag coefficient C as
D
a source term ∇.Σ in Eqn. 6). All the issues of principle obtainedusingIBMinthechannelgeometryforlargecylinder
arealreadycontainedinthesolutionofEqns.6and7for spacingL /L =4withFaxen’sanalyticalresultforflowpast
x y
purelyNewtonianflowinthegeometriesofinteresthere, a single cylinder in a channel [48]. Lower: Map of velocity
for which known benchmarks exist, as follows. magnitude |v| calculated using IBM (left) and from Ref. [49]
SanganiandAcrivos[45]derivedanalyticalexpressions (right) for a cylinder radius R=0.25, cylinder spacing Lx =
0.625andthroughputrateQ=1.0. Thecolourschemeofthis
for the drag coefficient C as a function of the cylinder
D flowmapisdifferentfromthatoftheothersinthismanuscript
radius R in Newtonian flow past the biperiodic array of
to match that of the original work [49] against which we are
Fig. 1 (upper), separately for the small cylinder limit
benchmarking.
R→0, and for the limit R→L /2=1/2 (in our units)
y
in which adjacent cylinders approach contact. These are
shown bythe dottedand dashedlines respectively inthe
upper panel of Fig. 3. Numerical results obtained us- of cylinder radius R for a fixed (wide) horizontal cylin-
ing our IBM method are in excellent agreement with derseparationLx =4.0inthechannelgeometryofFig.1
these predictions, as shown by the circles in the same (lower)areinexcellentagreementwiththeanalyticalpre-
panel. Also shown (by the diamonds, squares and tri- diction of Faxen [48] for flow past a single cylinder in a
angles) are results for the drag coefficient additionally channel. Fig. 1 (lower) shows that our results for the
obtained in the present study using two different numer- full flow field in the channel geometry are in excellent
ical methods that are independent of the IBM described agreement with those of Ref. [49].
above: a phase field method and a circular-propagator
method [46]. As can be seen, these numerical results for
the drag coefficient are in excellent agreement between
all our three methods across the full range of cylinder
VI. RESULTS: VISCOELASTIC FLOW PAST A
radii. Foroneparticularvalueofthecylinderradius, the
BIPERIODIC ARRAY OF CYLINDERS
fullflowfieldascalculatedinabiperiodicarrayisshown
in the lower panel of Fig. 3: on the left using the IBM,
andontherightusingthepropagatormethod. Excellent Having carefully benchmarked our codes against
agreement is seen between the two methods. known results in the literature, we now present our new
In Fig. 4 (upper), we show that numerical results ob- results for viscoelastic flow past the biperiodic array of
tainedusingourIBMforthedragcoefficientasafunction cylinders as sketched in Fig. 1 (top).
8
20
r = R + 0.05
)15
5
0
0.
+
R 10
=
r
(
D
FIG. 5. A schematic representation of the three flow types: λ 5
left: rotation, middle: simple shear, and right: extension
(which is also called pure shear).
0
0 0.1 0.2 R 0.3 0.4
A. Character of the flow field
FIG. 6. Scaling with cylinder radius R of the rate of defor-
mationλ inaNewtonianfluid(redsolidline)intheregions
D
In Fig. 2 above we presented the stationary constitu-
just above and below the cylinder where the nature of the
tivecurvesoftheOldroydB,fene-Pandfene-CRmodels, flow is simple shear. The plotted values of λ are taken at a
D
separately for the idealised flow fields of homogeneous distancer=R+0.05verticallyfromeachcylindercentre(see
simple shear and homogeneous planar extension. A key red points in schematic in top left). The black dash-dotted
aim of this work is to examine whether the thickening line shows the scaling 2V/(1−2R)2.
observed as a function of increasing Weissenberg num-
ber in the flow of a viscoelastic fluid through a porous
medium can be understood in terms of the fluid’s under-
lyingconstitutivebehaviourinthosesimplerprotocolsof r = R + 0.025
4
homogeneous simple shear and planar extension. 5)
2
In any porous flow geometry, however, the flow field 0
will of course vary in space (and sometimes also in time) 0.
+
in a complicated way, and will in general comprise an R
admixture of both shear and extensional components at = 2
any location. To quantify the flow field at any location, (r D arepcpirecaurlaanticneg of
therefore,wewritethelocalvelocitygradienttensorK as λ zones
the sum of symmetric and antisymmetric contributions,
D = 1(K+KT)andΩ= 1(K−KT),witheigenvalues
λ =2(cid:113)1D :D and λ =2(cid:113)1Ω:Ω respectively. Here 00 0.1 0.2 R 0.3 0.4
D 2 Ω 2
λ measurestherateofdeformation,andaccordinglythe
D FIG.7. ScalingwithcylinderradiusRoftheratedeformation
strengthofthedeformingeffectthattheflowisexpected
λ in a Newtonian fluid (red solid line) in the regions just
D
tohaveonthepolymerchains. λ characterisestherate
Ω fore and aft the cylinder where the nature of the flow is pure
of rotation, which does not deform the polymers. shear (i.e., extension). The plotted values of λ are taken
D
Out of these two eigenvalues we also construct the at a distance r = R+0.025 horizontally from each cylinder
frame-invariant, rate-independent parameter [6, 22] centre (see red points in schematic in top right). The black
dash-dotted line shows the scaling 2V/(πR).
λ2 −λ2 3
q = D Ω. (19)
λ2 +λ2
D Ω
This quantifies the nature of the flow field at any loca- studies, as described in Sec. V above. In examining our
tion, in the following way. A value q = +1 corresponds numerical results, however, we find that this change is
to pure extensional flow, which is sometimes also called relatively modest. Because of this, we shall frame the
pure shear flow (as expressed relative to Cartesian axes discussion of our results in trying to understand how the
inEqn14). Forq =−1theflowispurelyrotational,and polymer responds, with increasing Wi, to a velocity field
will have no deforming effect on the polymer chains. For that is assumed to be unchanged from that at Wi = 0.
q =0therateofstrainingisequaltotherateofrotation, Accordingly, we focus the discussion in this section on
givinganequalsuperpositionofrotationandpureshear. the flow field that pertains at Wi=0.
This corresponds to simple shear flow (as expressed rel- As can be seen in the left two columns of Fig. 9, the
ative to Cartesian axes in Eqn. 13). These three cases rateofdeformationλ isstrongestroundtheupperand
D
q =−1,0,+1 are sketched in Fig. 5. lower edges of the cylinder for all values of the cylinder
Inpractice,ofcourse,theflowfieldinanygivengeom- radiusR,andisdominatedbysimpleshear. Asshownin
etrywillchangeasafunctionofWeissenbergnumberWi. Fig. 6, the deformation rate in these regions scales with
Indeed, this effect is fully accounted for in our numerical cylinder radius R as VL/(L−2R)2. Based on this, we
9
we define a Weissenberg number
8
r = R + d
Vτ
Wi = . (22)
) 2 (L−2R)
d 6
+
R In any regime where the pressure drop is dominated by
=
r 4 thesqueezingofthefluidthroughthegapbetweenverti-
(
D callyadjacentcylinders,weexpectWi tobetherelevant
λ d = 0.05 2
Weissenberg number to characterise the flow. Note that
2 d = 0.1
d = 0.15 whilethescalingsprovidedbyFigs.7-8aretakenatrep-
resentative points in the fluid, they can only expected
0 to hold to within an order 1 prefactor. However we will
0 0.1 0.2 R 0.3 0.4
show that each of the Weissenberg number definitions
Eqs. 21,22 accurately captures the onset of thickening in
FIG. 8. Scaling with cylinder radius of the rate of deforma-
the regime of R for which it is expected to apply, justi-
tionλ inaNewtonianfluidintheregionscentredalongthe
D fying our choices.
diagonal45◦,i.e.,asthefluidjuststartstosqueezeinto(and
Having discussed the nature of the flow field in the
subsequentlymoveoutof)thecontractionbetweenvertically
adjacentcylinders. Theplottedvaluesofλ aretakenatdis- limit of Newtonian flow, Wi → 0, we now proceed to
D
tances r =R+d along any diagonal angled at 45◦ from the describe our numerical results for the flow response of
cylinder centre (see red points in schematic in top left). All the Oldroyd B, fene-P and fene-CR models as the Weis-
fourpointsareequivalentduetothesymmetryoftheNewto- senberg number Wi (according, in any regime, to the
nian flow field. The black dash-dotted line shows the scaling most relevant of the above definitions in that regime)
2.5V/(1−2R). increases with the polymer relaxation time τ.
define a Weissenberg number
B. Oldroyd B model
VLτ
Wi = . (20)
0 (L−2R)2
InFig.9weexploretheflowresponseoftheOldroydB
If any regime exists in which the pressure drop is dom- model as a function of the cylinder radius R (down each
inated by these regions of simple shear just above and column) and the adimensional polymer relaxation time
below the cylinder, we would expect Wi to be the rele- τ, proportional to the Weissenberg number (along each
0
vant Weissenberg number to characterise that regime. row). (In each case the cell height Ly = Lx ≡ L = 1.0
For the low porosity geometries with small R, there andthevelocitycharacterisingthethroughputV =1.0.)
alsoexistsaregionofreasonablystrongdeformationfore Shown in each column (beyond the first two) of Fig. 9
andaftofthecylinder,whichisextensionallydominated. is a colourmap of the largest eigenvalue λW of the poly-
To characterise this, we define a Weissenberg number mer conformation tensor W, which characterises the de-
gree to which the polymer molecules are deformed by
Vτ
Wi = . (21) the flow. (W = δ in a fluid equilibrated at rest.) For
1 πR small Weissenberg number (third and fourth columns),
This takes as its characteristic time the residence time the colourmap of the polymer deformation λ is essen-
W
of the polymer near the cylinder, πR/V. The cylinders tiallythesameasthatofthedeformationrateoftheflow
are widely spaced in this regime, so we expect the inter- field λ (first column), consistent with the fact that the
D
cylinder spacing to be a less important lengthscale in solution of the Oldroyd B model in the Newtonian limit
comparison. AsshowninFig.7,thisdefinitionprovidesa τ →0 is W =δ+2τD.
reasonableapproximationofthedeformationrateinthis Moving rightwards across the montage from the third
region. Therefore in any regime in which the pressure column, τ (and so Wi) progressively increases. For each
drop is dominated by these regions of extension fore or cylinderradiusR(alongeachrow),thefirstnoticeableef-
aftofthecylinder,weexpectWi tobetherelevantWeis- fect with increasing τ is that the bright regions of strong
1
senberg number to characterise the flow in that regime. polymer deformation associated with the regions of sim-
NotethatwhileWi shouldstrictlybelabelledasaDeb- ple shear along the top and bottom of the cylinder shift
1
orah number [51], for simplicity here we retain the label slightly downstream (rightwards). We suggest that this
Wi . arisesduetotheonsetwithτ oftheinfluenceofthecon-
1
Finally, there exists a region of moderately strong ex- vectivetermv.∇W inthepolymerconstitutiveequation.
tensional flow centred around the 45◦ diagonal lines, as- (In contrast, for τ =0 the convection term plays no role
sociated with the fluid just starting to squeeze into, and andthedynamicsarepurelylocal.) Thepolymerconfor-
then subsequently move out of, the gap between verti- mation at any location is then affected not only by the
callyadjacentcylinders. AsshowninFig.8,thedeforma- flow field at that location, but also receives information
tionrateinthisregionscalesasV/(L−2R). Accordingly, about the flow field immediately upstream (leftwards).
10
/ max
W
+1 +1
_ 102
max
W
101
0 -1
100
FIG. 9. Montage of flow states for varying relaxation time and cylinder radius. Rows show cylinders with radii in the range
R = 0.05 → 0.35. The first column shows the maximum eigenvalue of the symmetrised velocity gradient tensor λ (which
D
we normalise to one), and the second column shows the flow character q (colourbars are shown bottom left). The remaining
columnsshowtheprincipaleigenvalueλ ofthepolymerconformationtensorW forτ =0.02,0.05,0.1,0.2,0.4,0.8,1.0,1.6on
W
a log colourscale (colourbar in bottom right).
Corresponding to the montage of Fig. 9, curves of the tensor, focused in a wake along the centreline aft of the
normalised drag coefficient χ as a function of increas- cylinder. This is due to the region of extensional flow to
ing Wi are shown in Fig. 10, separately for each value the right of the stagnation point at the centre-aft point
of R (each row of the montage). The signature in the of the cylinder edge, which causes a strong extensional
dragcoefficientoftheinitialslightshiftrightwardsofthe stretching of the polymer molecules as the Weissenberg
flowpatternjustdecribedisaninitialdecreaseofχwith numberincreases. (Recallthisregionofextensionalflow,
Wi, relative to the value χ=1.0 in the Newtonian limit q = 1, in the first two columns of the montage for small
Wi→0. Intheexistingliteraturethiseffectissometime R. As noted above, the flow map itself also changes
suggested to stem from a shear thinning effect [22, 52]. with Wi, but to a relatively modest extent.) This exten-
However the Oldroyd B model does not shear thin, so sionalstretchingofthepolymerchainsmanifestsitselfas
that explanation cannot be valid here. We suggest in- a strong upturn in the drag coefficient following a mini-
stead that this initial decrease in drag arises due to the mumatWiup. AscanbeseeninFig.10a),thisoccursat
effect of the convective term in slightly (spatially) ‘de- a(nearly)constantvalueoftheWeissenbergnumberWi
1
laying’thebuild-upofviscoelasticstressatanylocation, for several different (small) values of the cylinder radius
relativetotheNewtoniancase. Inthisway,theregionsof R. This confirms that Wi , which we recall is defined by
1
strongest polymer deformation are shifted slightly away considering the rate of extensional deformation aft (and
from the regions of strongest shear rate. fore) of each cylinder, is indeed the relevant dimension-
MovingfurtherrightacrossthemontageofFig.9,more less rate to characterise the flow in this regime of low
dramatic changes in the flow field are evident. For small cylinder radius and high medium porosity.
values of the cylinder radius R, as τ (and so Wi) in- Jamesetal. studiedexperimentallytheflowofaBoger
creases rightwards across the montage, a bright streak fluid past a square array of cylinders in this regime of
develops in the colourmap of the polymer conformation small R/L = 0.09 → 0.18 [9]. As in our simulations,