Thickening of viscoelastic flow in a model porous medium

TL;DR

This study numerically investigates viscoelastic flow through a model porous medium, revealing mechanisms behind flow thickening as the Weissenberg number increases, with implications for understanding polymer solution flow behavior.

Contribution

It identifies two distinct mechanisms of flow thickening in porous media and explores how porosity influences these effects using numerical simulations.

Findings

Flow exhibits initial mild decrease then dramatic thickening with increasing Weissenberg number.

Two mechanisms of thickening are identified depending on obstacle spacing.

Flow remains steady across all Weissenberg numbers in a linear array configuration.

Abstract

We study numerically two-dimensional creeping viscoelastic flow past a biperiodic square array of cylinders within the Oldroyd B, FENE-CR and FENE-P constitutive models of dilute polymer solutions. Our results capture the initial mild decrease then dramatic upturn ('thickening') seen experimentally in the drag coefficient as a function of increasing Weissenberg number. By systematically varying the porosity of the flow geometry, we demonstrate two qualitatively different mechanisms underpinning this thickening effect: one that operates in the highly porous case of widely spaced obstacles, and another for more densely packed obstacles, with a crossover between these two mechanisms at intermediate porosities. We also briefly consider 2D creeping viscoelastic flow past a linear array of cylinders confined to a channel, where we find that the flow is steady for all Weissenberg numbers explored.

Equations

$$\frac{\Delta P}{L}=\frac{\eta V}{K},$$

$$C_D=\frac{\Delta P{L}_y}{\eta V},$$

$$\chi =\frac{C_D(\mathrm{Wi})}{C_D(\mathrm{Wi}\to 0)}.$$

$$\mathbf{T}=\mathbf{\Sigma }+2\eta \mathbf{D}-p\mathbf{I}.$$

$$\nabla \cdot \mathbf{T}=0,$$

$$\nabla .\mathbf{\Sigma }+\eta {\nabla }^2\mathbf{v}-\nabla p=0.$$

$$\nabla \cdot \mathbf{v}=0.$$

$$\mathbf{\Sigma }=G\left[f(\mathbf{W})\mathbf{W}-g(\mathbf{W})\mathbf{I}\right],$$

$$\stackrel{\nabla }{W}=-\frac{1}{\tau }\left[f(\mathbf{W})\mathbf{W}-g(\mathbf{W})\mathbf{I}\right]+\frac{{\ell }^2}{\tau }{\nabla }^2\mathbf{W},$$

$$\stackrel{\nabla }{W}\equiv \left({\partial }_t+\mathbf{v}\cdot \nabla \right)\mathbf{W}-{\nabla \mathbf{v}}^T\cdot \mathbf{W}-\mathbf{W}\cdot \nabla \mathbf{v}$$

$$f(\mathbf{W})=g(\mathbf{W})=\alpha (\mathbf{W})\equiv \frac{1}{1-\frac{Tr(\mathbf{W})}{{\Lambda }^2}},$$

$$f(\mathbf{W})=\alpha (\mathbf{W})\mathrm{and}\mathrm{g}(\mathbf{W})=1.$$

$$\nabla \mathbf{v}{|}_{\mathrm{shear}}=\dot{\gamma }\left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right).$$

$$\nabla \mathbf{v}{|}_{\mathrm{ext}}=\dot{ϵ}\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right).$$

$$\mathbf{F}(\xi ,t)=-\kappa \left(\mathbf{X}(\xi ,t)-{\mathbf{X}}_0(\xi )\right),$$

$$\frac{\partial \mathbf{F}(\xi ,t)}{\partial t}=-\kappa \mathbf{U}(\xi ,t),$$

$$\mathbf{U}(\xi ,t)={\int }_{\Omega }\mathbf{v}(\mathbf{r},t)\delta (\mathbf{r}-\mathbf{X}(\xi ,t))d\mathbf{x}.$$

$$\mathbf{f}(\mathbf{r},t)={\int }_{\Omega }\mathbf{F}(\xi ,t)\delta (\mathbf{r}-\mathbf{X}(\xi ,t))d\xi .$$

$$8h{\delta }_P(r)$$

$$q=\frac{{\lambda }_{\mathbf{D}}^2-{\lambda }_{\mathbf{\Omega }}^2}{{\lambda }_{\mathbf{D}}^2+{\lambda }_{\mathbf{\Omega }}^2}.$$

$${\mathrm{Wi}}_0=\frac{VL\tau }{(L-2R{)}^2}.$$

$${\mathrm{Wi}}_1=\frac{V\tau }{\pi R}.$$

$${\mathrm{Wi}}_2=\frac{V\tau }{(L-2R)}.$$

$$\underset{\text{advection, }{\mathbf{T}}_{\text{a}}}{\underbrace{\left(\mathbf{v}\cdot \nabla \right)\mathbf{\Sigma }}}=\underset{\text{distortion, }{\mathbf{T}}_{\text{d}}}{\underbrace{{\nabla \mathbf{v}}^T\cdot \mathbf{\Sigma }+\mathbf{\Sigma }\cdot \nabla \mathbf{v}}}+\underset{\text{remaining, }{\mathbf{T}}_{\text{r}}}{\underbrace{2G\mathbf{D}-\frac{1}{\tau }\mathbf{\Sigma }}}.$$

$${\lambda }_{\mathrm{a}}=\sqrt{{\mathbf{T}}_a:{\mathbf{T}}_a},{\lambda }_{\mathrm{d}}=\sqrt{{\mathbf{T}}_d:{\mathbf{T}}_d},{\lambda }_{\mathrm{r}}=\sqrt{{\mathbf{T}}_r:{\mathbf{T}}_r},$$

$${\mathrm{Wi}}_{\max }=\max ({\mathrm{Wi}}_1,{\mathrm{Wi}}_2).$$

$${\int }_0^{L_x}\nabla .\mathbf{v}dx$$

$${\int }_0^{L_x}\left[{\partial }_x{v}_x(x,y)+{\partial }_y{v}_y(x,y)\right]dx$$

$$\left[v_x(L_x,y)-v_x(0,y)\right]+{\int }_0^{L_x}{\partial }_y{v}_y(x,y)dx$$

$${\partial }_y{\int }_0^{L_x}{v}_y(x,y)dx$$

$${\int }_0^{L_x}{v}_y(x,y)dx$$

Full text

I Introduction

Flows of polymeric melts and solutions exhibit a rich phenomenology that has been widely studied Larson2014 . Many such materials exhibit viscoelastic properties which dramatically affect their behaviour in processes used in their industrial application. Examples range from melt extrusion through coating process to flow in porous materials. The latter example has particular relevance in the oil and gas industry where viscoelastic solutions are routinely employed for enhanced oil recovery, matrix stimulation and fracturing. For enhanced oil recovery, aqueous polymer solutions have long been used to control the Saffman-Taylor instability Homsy1987 ; Bensimon1986 that results from the displacement of viscous oil by lower viscosity brine within a porous reservoir rock. Since the earliest experiments it has been known that for certain polymeric solutions, even for single-phase flow, an anomalously high pressure gradient is observed when compared with the flow of an appropriate equivalent Newtonian solution. Whereas this phenomenon has been long known, no clear understanding of the root cause of the high observed drag has been forthcoming. In particular there has been no clear elucidation of the relative contributions of shear, extensional and elastic stresses. The excess pressure gradient directly impacts the industrial process by limiting the rate of injection of polymeric solution, where an excessive pressure gradient will lead to fracturing of the rock in the vicinity of the injector and to degradation of the polymer.

Experimentally, several authors have studied the flow of a viscoelastic fluid past a biperiodic array of obstacles arranged in two spatial dimensions Chmielewski1992 ; Skartsis1992 ; Khomami1997 ; Howe2015 ; James2012 ; James2016 . (The upper panel of Fig. 1 shows a sketch of the simplified cartoon of such a geometry that we shall study numerically in this work.) The dominant physical effect consistently reported in these experiments is a dramatic upturn in the adimensional pressure drop, known as the drag coefficient, relative to that for a Newtonian fluid of matched viscosity, for Weissenberg numbers $\mathrm{Wi}$ exceeding a critical value $\mathrm{Wi}_{\textrm{c}}$. (The Weissenberg number is the product of a characteristic shear-rate in the flow and polymer relaxation time $\tau$. Later in the text we define three Weissenberg numbers suited to the problem at hand.) This upturn is often also accompanied by the development of time-dependent flow fields, with crossing streaklines and structure in the third spatial dimension Chmielewski1992 ; Khomami1997 , into the page in the simplified sketch of Fig. 1. For low Weissenberg numbers a subtler effect is sometimes also seen, in which the pressure drop initially decreases slightly relative to the Newtonian case Khomami1997 before the upturn just described for $\mathrm{Wi}>\mathrm{Wi}_{\textrm{c}}$.

Some of these studies Clarke2015 ; Clarke2015a correlate the observed upturn in the drag coefficient (described therein via an apparent viscosity) found for polymer solution flow in outcrop rock samples with flows of the same solutions observed in 2D microfluidic networks. These networks comprise “pores” and randomly sized “throats” on a grid rotated $45^{\circ}$ to the average flow direction. The experiments were subsequently extended Howe2015 to show how the onset of additional drag depended on solution parameters. In particular, the onset of thickening was found to be well characterized by a Weissenberg number derived using a characteristic apparent shear rate proportional to flow rate together with a molecular relaxation time (see Fig. 11 in that work). In these studies, the transition to a time-dependent state resembling turbulence is marked by the appearance of crossing streaklines, a 3D effect which isn’t possible in the current 2D study. Nuclear magnetic resonance studies Mitchell2016 demonstrated time-dependent flows within a full 3D pore network (a rock) via an effective diffusion constant measurement. The analysis developed in the following paper can also be applied to these geometries at 45 degree orientations. Our results (not discussed here) suggest that the cylinder size has less effect on the flow character $q$ (defined in Section VI.1) than for $0^{\circ}$ orientations. However we defer a full study in these rotated geometries to future work.

From a theoretical viewpoint, early attempts to understand flow in porous media adopted a coarse-grained approach, discarding microscopic details in favour of macroscopic properties. For Newtonian flow, Darcy darcy1856 proposed a relation between the pressure drop per unit length of material $\Delta P/L$ and the mean velocity scale $V$

[TABLE]

where $\eta$ is the fluid’s viscosity. The permeability $K$ is a constant that should depend only on the properties of the medium, but is unknown a priori. A relation between the permeability and the medium’s porosity $\epsilon$ (the ratio of free volume to total volume) was later proposed in the Blake-Kozeny-Carman equation Carman1937 , which proved successful in describing a range of simple flows. However the validity of this macroscopic approach remains largely limited to Newtonian flows.

In theoretically understanding viscoelastic creeping flows in porous media, much of the progress has been made computationally Talwar1992 ; Talwar1995 ; Souvaliotis1996 ; LiuThesis1997 ; Alcocer1999 ; Alcocer2002 ; Gillissen2013 . Early simulations of two-dimensional (2D) viscoelastic flow past a bi-periodic square Talwar1992 ; Talwar1995 or rectangular Souvaliotis1996 array of cylinders observed a slightly reduced pressure drop at low Weissenberg numbers compared to that of a Newtonian fluid of matched viscosity, as seen experimentally. However the dramatic upturn in the drag coefficient at high $\mathrm{Wi}>\mathrm{Wi}_{\textrm{c}}$, which is the dominant physical effect seen experimentally, was not captured in these early numerical works, presumably due to the restricted computational processing power available at the time. It was however later captured in simulations of the Oldroyd B and FENE-CR models, also in 2D biperiodic arrays LiuThesis1997 .

Alcocer et al. Alcocer1999 ; Alcocer2002 investigated 2D flow of the FENE-CR model past a biperiodic array of cylinders, with a particular focus on the dependence of the effective permeability on the cell aspect ratio, for a fixed area fraction of cylinders. They demonstrated a non-monotonic dependence of permeability on aspect ratio. For a fixed aspect ratio and area fraction, they reported in an initial increase in permeability with increasing $\mathrm{Wi}$, equivalent to the initial decrease in drag in other studies.

Gillissen Gillissen2013 simulated 2D flow of the FENE-P model past a biperiodic hexagonal array of cylinders. This study convincingly captured both the initial downturn and then significant upturn in the drag coefficient seen experimentally as a function of Weissenberg number. Analysing the flow field as a function of space in terms of regions of pure shear (which is the same as extension), simple shear and pure rotation, they demonstrated a predominance of shear regions at low $\mathrm{Wi}$, with a progressive increase in elongational regions with increasing $\mathrm{Wi}$. At high $\mathrm{Wi}$ the polymer conformation tensor was found to be fully extended, showing the importance of finite chain extensibility in this regime.

De *et al. *studied 3D flow of a FENE-P fluid past an array of cylinders (both with and without walls) De2016 . They found an elastic instability whereby recirculating regions in the cylinder wake break symmetry and form a 3D structure, which occurs at a Deborah number $De=\tau V/R$ consistent with Ref. Smith2003 . All their simulation runs attained a time-independent steady-state.

Besides the biperiodic geometries just discussed, significant efforts have also been devoted to understanding viscoelastic flow past a single cylinder or linear array of cylinders confined to a channel: experimentally McKinley1993 ; LiuThesis1997 ; Moss2010a ; Ribeiro2014 ; Zhao2016 , by linear stability analysis Smith2003 ; Sahin2008 , and by direct numerical simulation Ribeiro2014 ; Liu1998a ; Hulsen2005 ; Oliveira2005 ; Claus2013 ; LiuThesis1997 ; Vazquez-Quesada2012 ; Grilli2013 . The lower panel of Fig. 1 shows a sketch of the simplified cartoon of such a geometry that we shall study numerically in this work.

By studying the flow of Boger fluids past a single cylinder in a channel, McKinley *et al. *McKinley1993 observed a transition from steady 2D to steady 3D flow in the downstream wake. At higher flow rates, they found another instability where time-dependent velocity oscillations form in the wake region. Liu LiuThesis1997 considered flow past a single cylinder, and widely and closely spaced linear arrays of cylinders. For the single cylinder, they observed a mild downturn in the drag at moderate $\mathrm{Wi}$, followed by an upturn at larger $\mathrm{Wi}$ which was accompanied by a transition from steady 2D to steady 3D flow. In the linear arrays, the transition was from steady 2D to time-dependent 3D structure in both cases. Moss and Rothstein Moss2010a studied flow of wormlike-micelles (WLMs) past a single cylinder for several ratios of cylinder diameter to channel width. They observe a significant decrease in the normalised pressure drop as a function of $\mathrm{Wi}$ which was attributed to the shear-thinning properties of the fluids. Of the two fluids tested, only one exhibited an instability, which was attributed to breakdown of the WLMs in the extensional flow in the cylinder wake. Using flow-induced birefringence measurements, the authors showed that shear flows at the channel walls were not necessary to produce the wake instability. Recent experiments on WLMs Zhao2016 found that upstream vortices formed at much larger $\mathrm{Wi}\sim 10^{3}$, and unsteady flow downstream at even larger $\mathrm{Wi}\sim 10^{4}$.

Simulations of 2D creeping viscoelastic flow past a linear array of cylinders confined in a channel Liu1998a ; Hulsen2005 captured an initial mild decrease then subsequent upturn in the drag coefficient as a function of increasing Weissenberg number. All the states observed were however steady as a function of time. In 2D simulations of flow past a single cylinder confined in a channel, temporal oscillations in the size of a recirculating region that forms downstream of the cylinder are seen Oliveira2005 ; Claus2013 , with an associated slight increase in the drag coefficient compared with the time-independent state.

Attempts to understand the onset at high $\mathrm{Wi}$ of the 3D time-dependent states seen experimentally for creeping viscoelastic flow past an array of cylinders in a channel have been made by performing a linear stability analysis for the dynamics of small amplitude 3D perturbations to an initially 2D flow state. By such an analysis, Smith et al. Smith2003 reproduced some of the experimentally observed characteristics of the instability, in particular the wavevector of the most unstable mode and the critical $\mathrm{Wi}$ at which the instability first arises. Sahin and Wilson Sahin2008 demonstrated that the wavelength of the instability scales with cylinder spacing for closely spaced cylinders, and with the size of the wake behind the cylinder for wider cylinder spacing. Both studies warn that this 3D instability potentially restricts the range of $\mathrm{Wi}$ over which purely 2D simulations might remain valid. Also we stress that linear stability does not preclude the existence of nonlinear instabilities, e.g., as discussed by Pan *et al. *Pan2013 .

Vázquez-Quesada and Ellero performed 2D smoothed particle hydrodynamics (SPH) simulations of viscoelastic flow past an array of cylinders in a channel Vazquez-Quesada2012 ; Grilli2013 , capturing the initial downturn and subsequent upturn in the drag coefficient as a function of Weissenberg number, as seen experimentally. At higher $\mathrm{Wi}$ the results of this study departed from other 2D numerical works LiuThesis1997 in reporting a transition to a time-dependent state, which the authors interpreted as viscoelastic turbulence. We do not find this viscoelastic turbulence, and discuss carefully the differences between our study and Refs. Vazquez-Quesada2012 ; Grilli2013 that might potentially explain this apparent discrepancy between the two studies.

Ribeiro et al. Ribeiro2014 studied fully 3D viscoelastic flow past a cylinder confined to a channel, both experimentally and numerically, for both a shear thinning and a Boger fluid. For the shear-thinning fluid they reported an elastic instability setting in upstream of the cylinder at critical value of $\mathrm{Wi}$ that depends on the cylinder height in the vorticity direction. For Weissenberg numbers just beyond onset of the instability, the system’s state is asymmetric and time-independent. A subsequent transition to a time-dependent state was reported at larger $\mathrm{Wi}$ still.

Among the simulation studies of viscoelastic flow just surveyed Talwar1992 ; Talwar1995 ; Souvaliotis1996 ; LiuThesis1997 ; Alcocer1999 ; Alcocer2002 ; Gillissen2013 , each considered one (or in some cases two) fixed ratio(s) of obstacle size to obstacle spacing, i.e., a fixed medium porosity. A key contribution of this work is systematically to vary the medium porosity over a broad range, from the limit of widely spaced obstacles to ones that nearly touch. In doing so, we shall demonstrate two qualitatively different mechanisms underpinning the dramatic upturn of the drag coefficient with Weissenberg number seen experimentally: one at low obstacle area fraction and another at high area fraction, with a crossover between these mechanisms at intermediate area fraction.

In neither case, however, do we find the thickening to be associated with the onset of a time-dependent flow, as often reported experimentally. The same is true in the 2D biperiodic simulations in the earlier literature Talwar1992 ; Talwar1995 ; Souvaliotis1996 ; LiuThesis1997 ; Alcocer1999 ; Alcocer2002 ; Gillissen2013 . One possible explanation for this discrepancy is that the time-dependent state seen experimentally is 3D in nature, and cannot be captured in purely 2D simulations.

The paper is structured as follows. In Secs. II and III we introduce the flow geometries and constitutive models to be studied. The governing parameters and dimensionless groups are summarised in Sec. IV. In Sec. V we outline our numerical methods, and provide benchmarks to validate them against known results in the Newtonian limit. We then present our results for viscoelastic flow: in Sec. VI for a biperiodic array of cylinders and in Sec. VII for a linear array of cylinders in a channel.

II Flow geometries

We shall study two different flow geometries. The first comprises a biperiodic array of cylinders, sketched in the upper panel of Fig. 1. The second comprises a linear array of cylinders in a channel bounded by solid walls, sketched in the lower panel of the same figure. In each case the flow cell has length $L_{x}$ horizontally and height $L_{y}$ vertically, and the cylindrical obstacle has radius $R$. At the boundaries of the cells represented by dashed lines all flow variables are ascribed periodic boundary conditions. At the cell boundaries represented by solid lines, and at the cylinder surface, conditions of no-slip and no-permeation apply. In each case we assume the flow to be translationally invariant into the page in Fig. 1, and simulate the flow in the two dimensions of the page only.

The flow, which we assume to be from left to right, can be imposed in two different ways. In the first, a given throughput $Q$ per unit time is prescribed (such that the characteristic velocity scale is then $V=Q/L_{y}$), with the pressure drop $\Delta P$ measured in response. Alternatively, we can impose the pressure drop $\Delta P$ and measure the resulting throughput $Q$.

In either case, a key experimental observable is the drag coefficient

[TABLE]

which measures the pressure drop normalised by the characteristic velocity scale and the solvent viscosity. In the limit of Newtonian flow this quantity depends only on the flow geometry. In non-Newtonian flow it also depends via the Weissenberg number on the nonlinear constitutive behaviour of the fluid in question. In reporting our numerical results below we shall typically show the drag coefficient at any given Weissenberg number, $C_{D}(\mathrm{Wi})$, normalised by the corresponding value in the Newtonian limit of zero Weissenberg number, $C_{D}(\mathrm{Wi}\to 0)$, defining

[TABLE]

III Constitutive models

We write the total stress $\bm{{T}}(\bm{{r}},t)$ in a fluid element at position $\bm{{r}}$ and time $t$ as the sum of a viscoelastic contribution $\bm{{\Sigma}}(\bm{{r}},t)$ from the polymer chains, a Newtonian solvent contribution of viscosity $\eta$, and an isotropic contribution with a pressure $p(\bm{{r}},t)$:

[TABLE]

The symmetric strain rate tensor $\bm{{D}}=\frac{1}{2}(\bm{{\nabla v}}+\bm{{\nabla v}}^{T})$ where $\bm{{\nabla v}}|_{\alpha\beta}=\partial_{\alpha}v_{\beta}$ and $\bm{{v}}(\bm{{r}},t)$ is the fluid velocity field.

Throughout we consider the creeping flow limit of zero Reynolds number. Here the condition of force balance requires the stress field $\bm{{T}}(\bm{{r}},t)$ to be divergence free:

[TABLE]

such that

[TABLE]

The pressure field $p(\bm{{r}},t)$ is determined by enforcing flow incompressibility:

[TABLE]

The dynamics of the polymeric stress $\bm{{\Sigma}}$ is specified by a viscoelastic constitutive model. In this work we consider three different phenomenological constitutive equations: the Oldroyd B, FENE-CR and FENE-P models Larson1988 . These each describe a dilute polymer solution by representing each polymer chain as a simplified dumbbell comprising two beads connected by a spring. The conformation tensor $\bm{{W}}=\langle\bm{{R}}\bm{{R}}\rangle$ is defined as the ensemble average $\langle\rangle$ of the outer dyad of the dumbbell end-to-end vector $\bm{{R}}$, which is taken to have unit length in equilibrium. The conformation of the polymer chains determines the viscoelastic stress according to

[TABLE]

with a constant modulus $G$. In addition to the spring force, each bead also experiences viscous drag against the solvent Larson1988 and stochastic thermal fluctuations. The conformation tensor is then taken to obey

[TABLE]

with a characteristic relaxation time $\tau$, where

[TABLE]

is the upper convected derivative Larson1988 . We have included a modification to the original equations by introducing a diffusive term, where $\ell$ is a small lengthscale below which gradients in $\bm{{W}}$ are attenuated. Similar modifications have been made to the Johnson-Segalman model in the context of shear banding Johnson1977 ; Radulescu1999 . Specifically in the context of porous media, Gillissen included diffusive terms in his study of a FENE fluid past a biperiodic array of cylinders Gillissen2013 . 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 Thomases2007 ; Thomases2011 . Without diffusion included, we find that the severe space- and time-step requirements limit the range of Weissenberg numbers that can feasibly be explored.

Gradient terms in $\bm{{W}}$ require a boundary condition at the walls and at the cylinder surface. For simplicity we choose zero gradient at the walls $\partial_{y}\bm{{W}}=0$ (when present). While the zero gradient condition can also explicitly be imposed on the cylinder surface111The method is analogous to the procedure used to enforce no-slip boundary conditions on the cylinder surface, see the discussion preceding Eq. 15. Essentially, we add a term to restore the gradient normal to the cylinder to zero., we find that our numerical method (which solves for $\bm{W}$ everywhere in the computational domain, both inside and outside of the cylinder) in fact naturally produces an emergent zero gradient at a typical radius $R+\ell$. We therefore find that it is sufficient to ensure that $\ell$ is small compared to any other physical lengthscale in the problem: we fix $\ell=0.01$ in all that follows and we have checked that our results are qualitatively (and in almost all cases quantitatively222The only exception is in the closely spaced cylinders case shown in Fig. 17 (right) where the drag upturn is slightly less steep for $\ell=0.005,0.0025$. The point of upturn in the drag is unchanged in all cases, and all states remain time-independent. ) unchanged by further decreasing $\ell$.

The functions $f(\bm{{W}})$ and $g(\bm{{W}})$ in Eqns. 8 and 9 are different in the three different models. The Oldroyd B model has $f(\bm{{W}})=g(\bm{{W}})=1$, which corresponds to assuming that the spring of each dumbbell is Hookean. For a sustained imposed extensional strain rate $\dot{\epsilon}>1/2\tau$ this model displays an extensional catastrophe in which the dumbbells stretch out indefinitely and the extensional stress grows indefinitely.

The phenomenological FENE-CR and FENE-P models regularise this catastrophe by insisting that the extensional stress of the polymer chains (dumbbells) must remain finite at all deformation rates. They do so by replacing the Hookean spring law by a non-linear law with finite-extensibility Chilcott1988 . The FENE-CR model has

[TABLE]

while the FENE-P model has

[TABLE]

The parameter $\Lambda$ in Eqn. 11 characterises the maximum extent to which any dumbbell can be stretched. We choose to express this in terms of $\delta=\Lambda^{-2}$, noting that the limit $\delta\to 0$ corresponds to Oldroyd B dynamics with infinite extensibility.

Under conditions of ideal viscometric simple shear flow, the imposed velocity gradient tensor in Cartesian ($x-y$) coordinates is

[TABLE]

In steady state under this applied flow, the shear stress $\Sigma_{xy}$ (which we also denote $\sigma$) as a function of $\dot{\gamma}$ defines the fluid’s shear constitutive curve, as shown in Fig. 2 for the three models considered here. The shear viscosity of the Oldroyd B and FENE-CR models is constant as a function of $\dot{\gamma}$. The FENE-P model shear thins to an extent determined by the value of $\delta$.

In an ideal viscometric planar extensional flow, the imposed velocity gradient tensor

[TABLE]

In steady state, the extensional stress $\Sigma_{xx}-\Sigma_{yy}$ (which we also denote $\sigma_{\rm E}$) as a function of $\dot{\epsilon}$ defines the fluid’s extensional constitutive curve, as shown in Fig. 2 for the three models. As can be seen, the extensional constitutive curve of the Oldroyd B model is undefined for flow rates $\dot{\epsilon}>1/2\tau$, consistent with our discussion above of the chain stretch catastrophe. In contrast, the FENE models have well-defined constitutive curves at all strain rates $\dot{\epsilon}$. Indeed for matched $\delta$ the two FENE models have the same extensional constitutive curves, with a finite limiting extensional viscosity $2G\tau/\delta$ as $\dot{\epsilon}\to\infty$.

IV Parameters and dimensionless groups

The eight parameters characterising the fluid, geometry, and imposed flow just described are: the solvent viscosity $\eta$, the polymer modulus $G$, the polymer viscoelastic relaxation timescale $\tau$, the polymer finite extensibility parameter $\delta$, the cell length $L_{x}$, the cell height $L_{y}$, the cylinder radius $R$ and the flow’s throughput rate $Q$.

We are free to choose units of mass, length and time, leaving five dimensionless groups as follows: the ratio of cylinder radius to gap height $\tilde{R}=R/L_{y}$, the ratio of cell length to cell height $\tilde{L}_{x}=L_{x}/L_{y}$, the ratio of solvent to total (zero shear) viscosity $\beta=\eta/\left(\eta+G\tau\right)$, the finite extensibility parameter $\delta$, and a Weissenberg number $\mathrm{Wi}$ characterising the strength of the velocity gradients compared to the inverse of the fluid’s stress relaxation time $\tau$. (We return below to define $\mathrm{Wi}$ precisely.) We drop tildes hereafter with the understanding that (for example) values for $R$ and $L_{x}$ are always quoted in units of $L_{y}$. To allow benchmarking of our results against the earlier literature we fix the viscosity ratio $\beta=0.59$ throughout. Remaining to be explored numerically are then the dimensionless cylinder radius $R$, the dimensionless cell length $L_{x}$, the finite extensibility parameter $\delta$ and the Weissenberg number $\mathrm{Wi}$.

The Weissenberg number is a dimensionless quantity characterising the scale of velocity gradients in the flow in units of the inverse relaxation time. On purely dimensional grounds, one simple possible definition is $V\tau/L_{y}$. However we have found it more useful to define two different Weissenberg numbers based on a characterisation of the velocity gradients that actually develop inside the flow geometry. We shall return to discuss these in Sec. VI below.

In our studies of the biperiodic geometry we consider a square cell with $L_{x}=L_{y}\equiv L=1.0$. For the channel geometry we take the cylinder spacing $L_{x}$ as an additional variable to be explored numerically.

V Numerical methods

We numerically solve the model equations 6, 7, 8 and 9 using a timestepping approach. Each main timestep comprises two separate substeps, as follows. In the first substep the viscoelastic stress $\bm{{\Sigma}}(\bm{{r}},t)$ is updated according to the viscoelastic constitutive equations 8 and 9 with a fixed velocity field $\bm{{v}}(\bm{{r}},t)$. This update is performed in the whole plane of Fig. 1, even inside the cylinder, where the velocity is zero to within the accuracy of our numerical approach. For this substep we adopt a method that we have used in our own previous work Fielding2006 , to which we refer the reader for details.

With that newly updated viscoelastic stress field, the 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 renders the procedure more complicated than in Ref. Fielding2006 (which considered a rectangular flow cell with no obstacles). In particular, the boundary conditions of no-slip and no-permeation must be satisfied for the fluid velocity round the edge of each obstacle. To tackle this we use an immersed boundary method (IBM)Peskin1972 ; Peskin2002 ; Mittal2005 , which couples a solution of Eqns. 6 and 7 on a regular Eulerian grid with an off-grid Lagrangian description of the cylinder surface. (We use the term Lagrangian to be consistent with the IBM literature, although in our particular problem the cylinder doesn’t advect with the flow.)

The cylinder surface is characterised by a one-dimensional curvilinear Lagrangian coordinate $\xi$ around it. (Although we use the word cylinder, in our 2D study it is of course represented by a circular cross section, with a 1D edge.) The location of a point at the Lagrangian coordinate $\xi$ is given in the Eulerian frame as $\bm{{X}}(\xi,t)$. A Lagrangian force density $\bm{{F}}(\xi,t)$ is then incorporated, calculated by prescribing the desired location of the cylinder $\bm{{X}}_{0}(\xi)$, and imposing a Hookean restoring force when the cylinder deviates from this:

[TABLE]

where $\kappa$ is a large spring constant. Given that the desired location $\bm{{X}}_{0}$ is independent of time, differentiating Eqn 15 with respect to time gives

[TABLE]

in which the Lagrangian velocity $\bm{{U}}$ at the cylinder surface is calculated from the velocity $\bm{{v}}(\bm{{r}},t)$ in the Eulerian frame as

[TABLE]

Eqn. 16 is evolved at each timestep using an explicit Euler algorithm. The Lagrangian force density $\bm{{F}}(\xi,t)$ then gives a force contribution in the Eulerian frame of

[TABLE]

This is incorporated as an additional source term to the left hand side of the generalised Stokes’ equation 6, which is then solved on a rectangular Eulerian grid in the full plane of Fig. 1 using the methods of Ref. Fielding2006 . Without walls, this can be done either imposing the overall throughput (and measuring the pressure drop) or imposing the pressure drop (and measuring the throughput), and we have checked that these are equivalent in all our biperiodic results presented below. With walls included, we always impose the latter quantity.

The transfer of information between the Lagrangian and Eulerian grids in Eqns. 17 and 18 is achieved in numerical practice by approximating the Dirac delta function $\delta(\bm{{r}})$ by a smoothed Peskin delta function $\delta_{P}(\bm{{r}})=\delta^{x}_{P}(x)\delta^{y}_{P}(y)$, in which  Lai2000

[TABLE]

Here $h=\Delta x=L_{x}/N_{x}=\Delta y=L_{y}/N_{y}$ is the Eulerian grid spacing, given a rectangular grid of $(N_{x},N_{y})$ points. The Lagrangian boundary of circumference $2\pi R$ is discretized into $M$ nodes of equal separation $\Delta s=2\pi R/M$. The optimal value of the ratio $\alpha=\Delta s/h$ is unknown a priori. Too small a value risks oversampling the boundary forces, while too large a value risks fluid leakage across the boundary. We set $\alpha=2$, and have checked that our results are robust to reasonable variations around this value. We have also ensured that all our results presented below are converged on the limit of grid spacing $h\to 0$ and of timestep $\Delta t\to 0$.

In the numerical solution just described, the methods employed to update the viscoelastic stress in the first substep of each timestep have been tested and benchmarked by ourselves in several previous publications Fielding2006 ; Fielding2007 . Therefore we focus here only on benchmarks to validate the second substep, in which the velocity field is calculated by solving Eqns. 6 and 7. The presence of viscoelasticity in this substep is trivial: it appears only as a source term $\nabla.\bm{{\Sigma}}$ in Eqn. 6). All the issues of principle are already contained in the solution of Eqns. 6 and 7 for purely Newtonian flow in the geometries of interest here, for which known benchmarks exist, as follows.

Sangani and Acrivos Sangani1982 derived analytical expressions for the drag coefficient $C_{D}$ as a function of the cylinder radius $R$ in Newtonian flow past the biperiodic array of Fig. 1 (upper), separately for the small cylinder limit $R\to 0$, and for the limit $R\to L_{y}/2=1/2$ (in our units) in which adjacent cylinders approach contact. These are shown by the dotted and dashed lines respectively in the upper panel of Fig. 3. Numerical results obtained using our IBM method are in excellent agreement with these predictions, as shown by the circles in the same panel. Also shown (by the diamonds, squares and triangles) are results for the drag coefficient additionally obtained in the present study using two different numerical methods that are independent of the IBM described above: a phase field method and a circular-propagator method HemingwayThesis2015 . 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 radii. For one particular value of the cylinder radius, the full flow field as calculated in a biperiodic array is shown in the lower panel of Fig. 3: on the left using the IBM, and on the right using the propagator method. Excellent agreement is seen between the two methods.

In Fig. 4 (upper), we show that numerical results obtained using our IBM for the drag coefficient as a function of cylinder radius $R$ for a fixed (wide) horizontal cylinder separation $L_{x}=4.0$ in the channel geometry of Fig. 1 (lower) are in excellent agreement with the analytical prediction of Faxen Faxens1946 for flow past a single cylinder in a channel. Fig. 4 (lower) shows that our results for the full flow field in the channel geometry are in excellent agreement with those of Ref. Ellero2011 . Finally, in Fig. 5 we compare profiles of the polymer stress $\Sigma_{xx}$ against previous work in the biperiodic geometry Souvaliotis1996 , again demonstrating excellent agreement.

VI Results: viscoelastic flow past a biperiodic array of cylinders

Having carefully benchmarked our codes against known results in the literature, we now present our new results for viscoelastic flow past the biperiodic array of cylinders as sketched in Fig. 1 (top).

VI.1 Character of the flow field

In Fig. 2 above we presented the stationary constitutive curves of the Oldroyd B, FENE-P and FENE-CR models, separately for the idealised flow fields of homogeneous simple shear and homogeneous planar extension. A key aim of this work is to examine whether the thickening observed as a function of increasing Weissenberg number in the flow of a viscoelastic fluid through a porous medium can be understood in terms of the fluid’s underlying constitutive behaviour in those simpler protocols of homogeneous simple shear and planar extension.

In any porous flow geometry, however, the flow field will of course vary in space (and sometimes also in time) in a complicated way, and will in general comprise an admixture of both shear and extensional components at any location. To quantify the flow field at any location, therefore, we write the local velocity gradient tensor $\bm{{\nabla v}}$ as the sum of symmetric and antisymmetric contributions, $\bm{{D}}=\tfrac{1}{2}(\bm{{\nabla v}}+\bm{{\nabla v}}^{T})$ and $\bm{{\Omega}}=\tfrac{1}{2}(\bm{{\nabla v}}-\bm{{\nabla v}}^{T})$, with eigenvalues $\lambda_{\bm{{D}}}=\sqrt{\tfrac{1}{2}\bm{{D}}:\bm{{D}}}$ and $\lambda_{\bm{{\Omega}}}=\sqrt{\tfrac{1}{2}\bm{{\Omega}}:\bm{{\Omega}}}$ respectively. Here $\lambda_{\bm{{D}}}$ measures the rate of deformation, and accordingly the strength of the deforming effect that the flow is expected to have on the polymer chains. $\lambda_{\bm{{\Omega}}}$ characterises the rate of rotation, which does not deform the polymers.

Out of these two eigenvalues we also construct the frame-invariant, rate-independent parameter Skartsis1992 ; Gillissen2013

[TABLE]

This quantifies the nature of the flow field at any location, in the following way. A value $q=+1$ corresponds to pure extensional flow, which is sometimes also called pure shear flow (as expressed relative to Cartesian axes in Eqn 14). For $q=-1$ the flow is purely rotational, and will have no deforming effect on the polymer chains. For $q=0$ the rate of straining is equal to the rate of rotation, giving an equal superposition of rotation and pure shear. This corresponds to simple shear flow (as expressed relative to Cartesian axes in Eqn. 13). These three cases $q=-1,0,+1$ are sketched in Fig. 6.

In practice, of course, the flow field in any given geometry will change as a function of Weissenberg number $\mathrm{Wi}$. Indeed, this effect is fully accounted for in our numerical studies, as described in Sec. V above. In examining our numerical results, however, we find that this change is relatively modest. Because of this, we shall frame the discussion of our results in trying to understand how the polymer responds, with increasing $\mathrm{Wi}$, to a velocity field that is assumed to be unchanged from that at $\mathrm{Wi}=0$. Accordingly, we focus the discussion in this section on the flow field that pertains at $\mathrm{Wi}=0$.

As can be seen in the left two columns of Fig. 7, the rate of deformation $\lambda_{\bm{{D}}}$ is strongest round the upper and lower edges of the cylinder for all values of the cylinder radius $R$, and is dominated by simple shear. As shown in Fig. 8, the deformation rate in these regions scales with cylinder radius $R$ as $VL/(L-2R)^{2}$. Because the flux is constant through any vertical slice, the mean velocity $V$ increases as it is forced through the narrow vertical gap between adjacent cylinders, producing an effective velocity that might be expected to scale as $V_{e}\sim VL/(L-2R)$. A typical shear-rate in the gap would then be $V_{e}/(L-2R)=VL/(L-2R)^{2}$. Based on this, we define a Weissenberg number

[TABLE]

If any regime exists in which the pressure drop is dominated by these regions of simple shear just above and below the cylinder, we would expect $\mathrm{Wi}_{0}$ to be the relevant Weissenberg number to characterise that regime.

For the low porosity geometries with small $R$, there also exists a region of reasonably strong deformation fore and aft of the cylinder, which is extensionally dominated. To characterise this, we define a Weissenberg number

[TABLE]

This takes as its characteristic time the residence time of the polymer near the cylinder, $\pi R/V$. The cylinders are widely spaced in this regime, so we expect the inter-cylinder spacing to be a less important lengthscale in comparison. As shown in Fig. 9, this definition provides a reasonable approximation of the deformation rate in this region. Therefore in any regime in which the pressure drop is dominated by these regions of extension fore or aft of the cylinder, we expect $\mathrm{Wi}_{1}$ to be the relevant Weissenberg number to characterise the flow in that regime. Note that while $\mathrm{Wi}_{1}$ should strictly be labelled as a Deborah number Colby2013 , for simplicity here we retain the label $\mathrm{Wi}_{1}$.

Finally, there exists a region of moderately strong extensional flow centred around the 45∘ diagonal lines, associated with the fluid just starting to squeeze into, and then subsequently move out of, the gap between vertically adjacent cylinders. As shown in Fig. 10, the deformation rate in this region scales as $V/(L-2R)$. Accordingly, we define a Weissenberg number

[TABLE]

In any regime where the pressure drop is dominated by the squeezing of the fluid through the gap between vertically adjacent cylinders, we expect $\mathrm{Wi}_{2}$ to be the relevant Weissenberg number to characterise the flow. Note that while the scalings provided by Figs. 9-10 are taken at representative points in the fluid, they can only expected to hold to within an order 1 prefactor. However we will show that each of the Weissenberg number definitions Eqs. 21,22 accurately captures the onset of thickening in the regime of $R$ for which it is expected to apply, justifying our choices.

Having discussed the nature of the flow field in the limit of Newtonian flow, $\mathrm{Wi}\to 0$, we now proceed to describe our numerical results for the flow response of the Oldroyd B, FENE-P and FENE-CR models as the Weissenberg number $\mathrm{Wi}$ (according, in any regime, to the most relevant of the above definitions in that regime) increases with the polymer relaxation time $\tau$.

VI.2 Oldroyd B model

In Fig. 7 we explore the flow response of the Oldroyd B model as a function of the cylinder radius $R$ (down each column) and the adimensional polymer relaxation time $\tau$, proportional to the Weissenberg number (along each row). (In each case the cell height $L_{y}=L_{x}\equiv L=1.0$ and the velocity characterising the throughput $V=1.0$.)

Shown in each column (beyond the first two) of Fig. 7 is a colourmap of the largest eigenvalue $\lambda_{\bm{{W}}}$ of the polymer conformation tensor $\bm{{W}}$, which characterises the degree to which the polymer molecules are deformed by the flow. ($\bm{{W}}=\bm{{\delta}}$ in a fluid equilibrated at rest.) For small Weissenberg number (third and fourth columns), the colourmap of the polymer deformation $\lambda_{\bm{{W}}}$ is essentially the same as that of the deformation rate of the flow field $\lambda_{\bm{{D}}}$ (first column), consistent with the fact that the solution of the Oldroyd B model in the Newtonian limit $\tau\to 0$ is $\bm{{W}}=\bm{{\delta}}+2\tau\bm{{D}}$.

Moving rightwards across the montage from the third column, $\tau$ (and so $\mathrm{Wi}$) progressively increases. For each cylinder radius $R$ (along each row), the first noticeable effect with increasing $\tau$ is that the bright regions of strong polymer deformation associated with the regions of simple shear along the top and bottom of the cylinder shift slightly downstream (rightwards). Given that this shift is in the direction of flow, one might expect that this arises due to the increasing influence of the advective term $\bm{{v}}.\nabla\bm{{W}}$ in the polymer constitutive equation with $\tau$. (In contrast, for $\tau=0$ the advection term plays no role and the dynamics are purely local.) The polymer conformation at any location would then be affected not only by the flow field at that location, but would also receive information about the flow field immediately upstream (leftwards).

However, elastic effects arising due to the distortion term cannot be neglected. In order to quantify the relative strength of advection and distortion, we first note that the constitutive equation Eq. 9 can be rewritten in terms of the polymer stress (as opposed to the conformation tensor $\bm{{W}}$). In steady state this reads

[TABLE]

The three tensorial contributions can be contracted to scalars as

[TABLE]

and integrated over space to yield a single value for each term. In Fig. 11 (top) we plot these integrals as a function of $\mathrm{Wi}_{2}$. This shows that for small to moderate values of $\mathrm{Wi}_{2}$, the steady-state equation is dominated by distortion rather than advection. This result is reinforced by the colourmaps in Fig. 11 (bottom) that show that distortion (near the cylinder tops) is stronger than advection (just before and after the narrowest vertical point). This suggests that the observed downturn in the drag is mainly a result of distortion rather than advection. Interestingly, beyond a certain value of $\mathrm{Wi}_{2}$, the advective term becomes comparable to distortion. We will later show that this coincides with the point at which the drag dramatically increases, $\mathrm{Wi}_{2}^{\rm up}$ (see black triangles in Fig. 12).

Corresponding to the montage of Fig. 7, curves of the normalised drag coefficient $\chi$ as a function of increasing $\mathrm{Wi}$ are shown in Fig. 12, separately for each value of $R$ (each row of the montage). The signature in the drag coefficient of the initial slight shift rightwards of the flow pattern just described is an initial decrease of $\chi$ with $\mathrm{Wi}$, relative to the value $\chi=1.0$ in the Newtonian limit $\mathrm{Wi}\to 0$. In the existing literature this effect is sometime suggested to stem from a shear thinning effect Moss2010 ; Gillissen2013 . However the Oldroyd B model does not shear thin, so that explanation cannot be valid here. We suggest instead that this initial decrease in drag arises due to the effect of the elastic distortion terms in slightly (spatially) ‘delaying’ the build-up of viscoelastic stress at any location (recall Fig. 11), relative to the Newtonian case. In this way, the regions of strongest polymer deformation are shifted slightly away from the regions of strongest shear rate.

Moving further right across the montage of Fig. 7, more dramatic changes in the flow field are evident. For small values of the cylinder radius $R$, as $\tau$ (and so $\mathrm{Wi}$) increases rightwards across the montage, a bright streak develops in the colourmap of the polymer conformation tensor, focused in a wake along the centreline aft of the cylinder. This is due to the region of extensional flow to the right of the stagnation point at the centre-aft point of the cylinder edge, which causes a strong extensional stretching of the polymer molecules as the Weissenberg number increases. (Recall this region of extensional flow, $q=1$, in the first two columns of the montage for small $R$. As noted above, the flow map itself also changes with $\mathrm{Wi}$, but to a relatively modest extent.) This extensional stretching of the polymer chains manifests itself as a strong upturn in the drag coefficient following a minimum at $\mathrm{Wi}^{\rm up}$. As can be seen in Fig. 12a), this occurs at a (nearly) constant value of the Weissenberg number $\mathrm{Wi}_{1}$ for several different (small) values of the cylinder radius $R$. This confirms that $\mathrm{Wi}_{1}$, which we recall is defined by considering the rate of extensional deformation aft (and fore) of each cylinder, is indeed the relevant dimensionless rate to characterise the flow in this regime of low cylinder radius and high medium porosity.

James studied experimentally the flow of a Boger fluid past a square array of cylinders in this regime of small $R/L=0.09\to 0.18$ James2012 . As in our simulations, they found the flow to be steady (up to values of their defined Deborah number $De\equiv\tau V/L\approx 4$). They likewise reported only a small downturn in the drag coefficient as a function of increasing $De$, before a pronounced upturn at $De=0.5$.

Returning to Fig. 7, for larger values of the cylinder radius $R$, the dominant effect as $\tau$ increases rightwards across the montage is that the layers of strong polymer deformation in the shear fields on the upper and lower edges of the cylinder intensify, and are supplemented by the development of secondary ‘layers’ of strong polymer deformation above and below the cylinder (one layer above the cylinder and one below it). These appear to originate in the contraction as the flow moves from the left hand edge of each snapshot into the narrow gaps between vertically adjacent cylinders, with these secondary ‘layers’ also advected downstream into the vertical gap. This again manifests itself as a strong upturn in the drag coefficient following a minimum at $\mathrm{Wi}^{\rm up}$. As seen in Fig. 12b), this upturn occurs at a constant value of $\mathrm{Wi}_{2}$ for several different values of the cylinder radius $R$. This confirms that $\mathrm{Wi}_{2}$ is indeed the relevant Weissenberg number to characterise the flow in this regime of larger cylinder radius and smaller medium porosity.

The values of $\mathrm{Wi}_{1}$ and $\mathrm{Wi}_{2}$ at the minima in the drag coefficient curves of Fig. 12 are plotted as a function of cylinder radius $R$ in Fig. 13. As can be seen, the upturn occurs at a roughly constant value of

[TABLE]

The crossover between which of $\mathrm{Wi}_{1}$ and $\mathrm{Wi}_{2}$ dominates at any value of $R$ occurs at $R_{\rm c}=L/(2+\pi)\approx 0.194$. For cylinder radii $R<R_{\rm c}$ we have $\mathrm{Wi}_{1}>\mathrm{Wi}_{2}$: in this regime, the effects of the extensional wake in the relatively wide horizontal gap between horizontally adjacent cylinders dominate those of shear in the gap between vertically adjacent cylinders. For $R>R_{\rm c}$ we have $\mathrm{Wi}_{2}>\mathrm{Wi}_{1}$: in this regime the squeezing of the fluid into the now narrower vertical gap between vertically adjacent cylinders dominates any effects in the now smaller horizontal extensional wake. This effectiveness of $\mathrm{Wi}_{\rm max}$ in characterising the flow across the full range of cylinder radii is also seen via the master plot of the drag coefficient $\chi$ as a function of $\mathrm{Wi}_{\rm max}$ in Fig. 12c): the upturn occurs at a fixed value of $\mathrm{Wi}_{\rm max}$ for all values of $R$.

VI.3 Fene models

We consider now the effect of finite dumbbell extensibility on the phenomena just discussed. In particular, we shall report the drag coefficient as a function of Weissenberg number in each of the FENE-P and FENE-CR models, for two different values of the cylinder radius: $R=0.1<R_{\rm c}\approx 0.194$ and $R=0.35>R_{\rm c}\approx 0.194$. The case $R=1/2.4\approx 0.41$ was considered previously LiuThesis1997 for both FENE models, and our findings for $R=0.35$ will be qualitatively consistent with that study.

Our results for the FENE-P model are shown in Fig. 14a). We recall from Fig. 2 that, under conditions of homogeneous viscometric flow, this model thins in both shear and extension. Consistently, we find that the drag coefficient is smaller for the FENE-P model ($\delta>0$) than for the Oldroyd B model ($\delta=0$) for both values of $R$ considered.

Our results for the FENE-CR model are shown in Fig. 14b). We recall from Fig. 2 that, under conditions of homogeneous viscometric flow, this model thins only in extension but not in shear. In the porous geometry studied here, for a cylinder radius $R=0.1<R_{\rm c}$ the drag coefficient is lower in FENE-CR ($\delta>0$) than in Oldroyd B ($\delta=0$). This is consistent with the extensional thinning of FENE-CR, and with the fact that the flow is dominated by the extensional wake aft the cylinder for this value of $R$. In contrast, for a cylinder radius $R=0.35>R_{\rm c}$, the drag coefficient is larger in FENE-CR than in Oldroyd B. Clearly, this observation lacks any obvious explanation in terms of the homogeneous constitutive curves of Fig. 2. Feasibly, it could arise because the finite dumbbell extensibility reduces the extent to which the molecules are stretched and reoriented as they transit the contraction flow en route into the gap between vertically adjacent cylinders, causing them then to confer a greater shear stress in that gap. We do not provide any evidence to support this claim, however.

VII Results: viscoelastic flow past an array of cylinders in a channel

We now present our results for the channel geometry sketched in Fig. 1 (lower), comprising a periodic linear array of cylinders bounded by solid walls. This has been studied widely in the existing literature Grilli2013 ; Liu1998a ; Smith2003 ; Oliveira2005 ; Sahin2008 ; Vazquez-Quesada2012 ; Sahin2013 . To allow a comparison between our results and some of those earlier studies we fix the cylinder radius $R/L_{y}=0.25$, as in Refs. Liu1998a ; Vazquez-Quesada2012 . This leaves the horizontal distance $L_{x}$ between cylinder centres as the geometrical parameter to be varied numerically. Noting that $L_{x}=2R=0.5L_{y}$ corresponds to touching cylinders and $L_{x}\to\infty$ to the limit of a single cylinder, we shall present results for two cases: closely spaced cylinders with $L_{x}/L_{y}=0.625$, and widely spaced cylinders with $L_{x}/L_{y}=1.5$.

In simulating this channel geometry we impose the pressure drop $\Delta P$ and measure the resulting flux $Q$. From these quantities, we can then calculate the normalised drag coefficient in Eqn. 3. In this walled geometry, the pressure drop equates to a sum comprising two contributions: one stemming from the drag on the cylinder plus another stemming from the drag on the wall. However, many earlier works in the literature report only the contribution from the drag on the cylinder. To allow a direct comparison with those works, we shall report a modified normalised drag coefficient $\bar{\chi}$, removing from Eqn. 2 the contribution to the pressure drop stemming from the force on the wall. We adopt from herein the literature definition $\mathrm{Wi}_{lit.}=\tau V/R$.

Flowmaps for the case of widely spaced cylinders with $L_{x}/L_{y}=1.5$ are shown in Fig. 15. As can be seen, the flow character as quantified by the parameter $q$ shows extensional regions fore and aft of the cylinder, regions of simple shear at the channel walls, and regions of simple shear just above and below the cylinder. As the Weissenberg number increases from left to right across the montage, the polymer conformation tensor becomes strongly deformed in the region of extensional wake aft of the cylinder. The corresponding normalised drag coefficient in Fig. 17 (left) shows similar behaviour as a function of Weissenberg number as in the biperiodic geometry, with an initial downturn then large upturn. This was reported also in the earlier studies of Liu1998a ; Vazquez-Quesada2012 .

Flowmaps for the case of closely spaced cylinders with $L_{x}/L_{y}=0.625$ are shown in Fig. 16. In this regime, the region between adjacent cylinders is effectively shielded from the main flow, as can be seen in the colourmap of $\lambda_{\bm{{D}}}/\lambda_{\bm{{D}}}^{max}$ and indeed the polymer conformation tensor remains largely undeformed in this region. The flow character as quantified by $q$ shows pronounced regions of extensional flow along the lines projecting diagonally outwards into the fluid from the cylinder centre, in the region where the fluid just starts to squeeze into (and subsequently move out of) the vertical gap between the cylinder and the channel walls. The corresponding drag coefficient (Fig. 17, right) shows a less pronounced downturn as a function of $\mathrm{Wi}$ than for the case of widely spaced cylinders, as also reported in Ref. Liu1998a .

After an initial startup transient, all of the states reported in Figs. 15, 16 and 17 reached a time-independent steady state. While this is consistent with most previous 2D studies LiuThesis1997 ; Liu1998a ; Hulsen2005 , it contradicts recent 2D works Vazquez-Quesada2012 ; Grilli2013 which reported time-dependent, turbulent-like states in 2D simulations of an Oldroyd-B fluid, in the same geometry as considered here. We now attempt to understand this discrepancy, first by appealing to the properties of incompressible 2D flow and then by discussing in turn the other differences between the two studies that might potentially explain this, including with regards to diffusivity, inertia, and fluctuations.

In Ref. Vazquez-Quesada2012 , the authors use the rms of the time-signal $\langle v_{y}(t)\rangle$ (where the average is taken over all space) as an order parameter for the transition to time-dependent states. Here we show that incompressible 2D flow with no-slip boundary conditions strictly requires $\langle v_{y}(t)\rangle=0$. This means that any numerical scheme based on incompressible hydrodynamics (such as that used in this paper) cannot hope to reproduce the fluctuations of Ref. Vazquez-Quesada2012 . Integrating the incompressibility condition along the length of the channel gives

[TABLE]

The left terms disappear due to periodic boundary conditions (alternatively because the flux cannot vary with $x$), leaving

[TABLE]

where $C$ is a constant independent of $y$. This integral must be zero at the boundaries (because $v_{y}=0$ at the walls) so $C=0$, meaning

[TABLE]

This can easily be generalised to include solid obstacles such as a cylinder, giving the same result. Therefore for incompressible hydrodynamics, $\langle v_{y}\rangle$ is a quantity that should be exactly zero, and not used as an order parameter. A possible explanation for the results of Ref. Vazquez-Quesada2012 could be that the flow in that case was slightly compressible, and the observed fluctuations (above Weissenberg numbers of order 1) were in the density.

The above analysis only pertains to fluctuations in the $y$-component of the velocity: fluctuations in the flux (and therefore $\mathrm{Wi}$) are permitted, as are time-dependent cylinder drag or lift forces. We now demonstrate that these allowable fluctuations only arise in our simulations if insufficient numerical resolution is used. (Recall that we shall also return below to discuss several other possible sources of the discrepancy between our work and that of Ref. Vazquez-Quesada2012 .) Focusing on the closely spaced geometry with $L_{x}/L_{y}=0.625$, which produced strongly fluctuating states in Ref. Vazquez-Quesada2012 , we repeat the simulations shown in Fig. 17 at reduced numerical resolution. In Fig. 18 we plot the standard deviation of the time series $\mathrm{Wi}(t)$ against the mean (after discarding the startup transient). For lowest resolution studied, we observe the onset of apparent time-dependent behaviour at $\mathrm{Wi}=\mathrm{Wi}_{c}\sim 1$. Similar to the results of Ref. Vazquez-Quesada2012 , the magnitude of the fluctuations approximately scales as $0.03(\mathrm{Wi}-\mathrm{Wi}_{c})^{1/2}$, suggestive of a Hopf bifurcation. For the larger resolution, $N_{y}=256$, we observe no fluctuations across the full range of $\mathrm{Wi}$ shown. However for the largest value of $\mathrm{Wi}$ shown in Fig. 17 (left), even $N_{y}=256$ is insufficient and similar fluctuations develop (not shown); these fluctuations disappear in our highest resolution simulations at $N_{y}=512$ and $N_{y}=768$.

Having shown that we do not find viscoelastic turbulence in our simulations, we return to discuss the several differences between out work and that of Refs. Grilli2013 ; Vazquez-Quesada2012 that might potentially explain this. First the model presented here includes a diffusive term. As discussed by Sureshkumar and Beris Sureshkumar1995 , this can have a stabilising effect, suppressing possible numerical instabilities. Here we have made an effort to minimise the effect of such a term by getting as close as is numerically possible to the dual limit $\Delta x,\Delta y/\ell\to 0$, $\ell/L_{x},L_{y}\to 0$, which ensures that the diffusive lengthscale $\ell$ is well resolved yet small compared to any features of the flow field. This approach suppresses spurious numerical instabilities whilst keeping physical ones (if present). A second explanation could be due to the presence of small but non-zero inertia in Ref. Vazquez-Quesada2012 . As discussed by Hoda *et al. *Hoda2008 , weak inertial effects can provide a mechanism by which elastic disturbances can be amplified. Finally it is conceivable that the time-dependent states observed in Ref. Vazquez-Quesada2012 are the result of a nonlinear instability. This mechanism has been proposed to explain observed instabilities in Poiseuille flow of viscoelastic fluids, which are believed to be linearly stable Bertola2003 ; Pan2013 . A final possibility concerns spatial resolution. While the number of SPH particles cannot necessarily be directly compared to the number of finite-difference grid points, we note that the second largest particle resolutions in Ref. Vazquez-Quesada2012 would roughly correspond to the lowest resolution grid in Fig. 18 for which we see fluctuations.

VIII Conclusions

In this work, we first studied two dimensional creeping flow of the Oldroyd B, FENE-CR and FENE-P viscoelastic fluids past a biperiodic square array of cylinders. Our aim has been to understand the dramatic upturn reported experimentally in the drag coefficient as a function of increasing Weissenberg number.

By performing simulations across a wide range of values of the porosity of the flow geometry, from ‘dilute’ to near-touching cylinders, we have demonstrated two qualitatively different mechanisms that may separately underpin this thickening effect. The first operates in the highly porous case of widely spaced cylinders, and involves a strong stretching of the polymer chains in the extensional wake aft of each cylinder. The second operates in the regime of more densely packed obstacles, and involves a strong deformation of the polymer chains as they squeeze into the vertical gap between vertically adjacent cylinders. Two different Weissenberg numbers separately characterise each of these regimes, and we have demonstrated that the upturn in the drag coefficient occurs at a fixed value of the maximum of these two numbers across the full range of medium porosities. We have also studied the creeping flow of an Oldroyd B fluid past a linear array of cylinders confined to a channel bounded by solid walls, where we have found that the flow remains steady for all Weissenberg numbers explored.

All the simulations in this work have assumed a purely two-dimensional flow, with translational invariance into the page in the sketches of Fig. 1. In experimental practice, the upturn in the drag coefficient is often accompanied by the onset of three dimensional flows (which are furthermore often time-dependent). To capture these effects, we would need to perform simulations in fully three dimensions. The two dimensional results reported here should clearly be treated with caution for Weissenberg numbers exceeding $O(1)$, where three dimensional effects may pertain.

James James2016 recently considered the ratio of shear- to extensionally-induced first normal stresses in viscoelastic flow past a biperiodic array of cylinders. They demonstrated an $O(1)$ lower bound on this, implying that shear generated first normal stresses $N_{1}$ cannot be neglected. Wagner and McKinley Wagner2016 recently examined the response of an Oldroyd-B fluid to a flow with a character-parameter $\alpha(t)$ varying sinusoidally in time, intended to mimic the variation experienced by a given fluid element as it moves through a periodic flow cell of the kind studied in this work. For large Deborah numbers $De$ (defined as the ratio of the polymer relaxation time to the time required to pass one unit cell), they demonstrated that attempts to predict the first normal stress $N_{1}$ just above and below the cylinder using the local Weissenberg number fails to reflect the previous flow history and drastically underestimates the value of $N_{1}$. In view of these recent works, in future numerical studies of the kind performed here it would clearly be interesting to consider more explicitly the role of normal stresses.

Acknowledgements

EJH and SMF gratefully acknowledge financial support for this work from Schlumberger Gould Research. We also thank the referees for their valuable comments on the paper.