Table Of Content

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 materials. 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ázquez-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 obstacles 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ázquez-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 dimensionless 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 finite 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 limiting extensional viscosity 2Gτ/δ as (cid:15)˙→∞. 2 10 /G σ IV. PARAMETERS AND DIMENSIONLESS e 0 GROUPS 10 The eight parameters characterising the fluid, geometry, 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 homogeneous 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 different 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 substep 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 obstacles). In particular, the boundary conditions of no-slip 0 and no-permeation must be satisfied for the fluid velocity 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 presented 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 calculated 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 number 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 degree 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,