Viscoelastic flow past mono- and bidisperse random arrays of cylinders: flow resistance, topology and normal stress distribution

We investigate creeping viscoelastic fluid flow through two-dimensional porous media consisting of random arrangements of monodisperse and bidisperse cylinders, using our finite volume-immersed boundary method introduced in S. De, et al., J. Non-Newtonian Fluid Mech., 2016, 232, 67-76. The viscoelastic fluid is modeled with a FENE-P model. The simulations show an increased flow resistance with increase in flow rate, even though the bulk response of the fluid to shear flow is shear thinning. We show that if the square root of the permeability is chosen as the characteristic length scale in the determination of the dimensionless Deborah number (De), then all flow resistance curves collapse to a single master curve, irrespective of the pore geometry. Our study reveals how viscoelastic stresses and flow topologies (rotation, shear and extension) are distributed through the porous media, and how they evolve with increasing De. We correlate the local viscoelastic first normal stress differences with the local flow topology and show that the largest normal stress differences are located in shear flow dominated regions and not in extensional flow dominated regions at higher viscoelasticity. The study shows that normal stress differences in shear flow regions may play a crucial role in the increase of flow resistance for viscoelastic flow through such porous media.


Introduction
Understanding the flow behavior of viscoelastic fluids through complex porous media is of prime importance due to its versatile applications, like enhanced oil recovery, catalytic polymerization, composite manufacturing, and in biological, geophysical, environmental and other processes. [1][2][3] Darcy's law 2,4 can be used to effectively characterize the creeping flow response of a Newtonian fluid to an imposed pressure drop or, conversely, the pressure drop associated with a certain imposed flow rate. However, when the fluid is non-Newtonian in nature, especially when it is viscoelastic, the flow behavior can no longer by predicted by Darcy's law and is often difficult to understand. 5 The difficulty is caused by the complex interactions and multiscale nature of the porous medium and fluid rheology. 6,7 There are ample experimental observations of an excess pressure drop, appearing beyond a critical flow rate of the viscoelastic fluid. 5,[8][9][10][11] This increased flow resistance has a strong correlation with the appearance of elastic instabilities in the viscoelastic fluid flow in the porous medium. 12,13 This may even apply on the scale of single pores: Chmielewski et al. 9 experimentally and numerically showed an increased pressure drop and departure from Darcy's law for flow of an elastic fluid through a single converging diverging geometry. Elastic instabilities and appearance of velocity and pressure fluctuations were reported in the experimental work of Arora et al. 12 The concept of elastic turbulence in microfluidic flow channels, and its relation to curvature of streamlines, has been reported extensively in recent decades. [14][15][16] More recently, the presence of elastic instabilities has also been confirmed in model porous media. [17][18][19] Most researchers have tried to link the excess pressure drop in these porous media to enhanced extensional flows and hysteresis in conformational stress. 20,21 However, our recent work on viscoelastic fluid flow through symmetric and asymmetric periodic sets of cylinders shows that most viscoelastic energy dissipation occurs in the shear-dominated regions of the flow domain, rather than extensional-flow-dominated regions. 22 Experimental findings of James et al. 23 also suggest that viscoelastic normal shear stresses are important in porous media flow. Thus, whether the increased flow resistance is mainly caused by extensional or shear flow mechanisms, is highly debated in the present literature.
Numerical studies can help in elucidating the dominant mechanism behind the enhanced pressure drop. Most numerical studies are limited to flow of power law fluids, which are inelastic non-Newtonian fluids. 24,25 The flow behavior of viscoelastic fluids through porous media is largely different from such power law fluids. Moreover, viscoelastic flow simulations have been mostly performed for regular linear arrays of cylinders or spheres. 26,27 Most of these numerical simulations were unable to fully capture the increased flow resistance, as compared to experimental findings on random porous media. [27][28][29][30] Liu et al. 27 employed a FENE based model to simulate viscoelastic fluid flow through an array of cylinders. However no such increase in flow thickening was observed. This might be due to the selection of a simpler model pore geometry, or due to limitations in the computational framework used in previous studies. Thus, the current literature shows that the complex interaction between the viscoelastic fluid rheology and pore structure is still not well understood. Moreover, numerical simulations of flow of viscoelastic fluids through more complex porous media are clearly missing.
In this paper, we report on numerical simulations of viscoelastic fluid flows through two dimensional disordered porous media consisting of packed arrangements of monodisperse and bidisperse sets of cylinders, using a finite volume and immersed boundary (FVM-IB) methodology. 1 We will compare in detail the flow structures and viscoelastic stresses in these two different sets of porous media, for different porosities and with increasing viscoelasticity (flow rate). We will show how the flow resistance changes with increasing Deborah number for all flow geometries. An analysis of flow topology enables us to understand how the shear, extension and rotational dominated flow regions are distributed in the porous media and evolve with increasing flow rate. Next a detailed analysis of the viscoelastic normal stress and its correlation to flow topology will give insight to a probable mechanism of the increased pressure drop and its relation to normal stress differences.

Constitutive equations
The fundamental equations for an isothermal incompressible viscoelastic flow are the equations of continuity and momentum, and a constitutive equation for the non-Newtonian stress components. The first two equations are as follows: Here u is the velocity vector, r is the fluid density (assumed to be constant) and p is the pressure. s is the viscoelastic polymer stress tensor. In the momentum equation, the Newtonian solvent stress is explicitly added and defined as 2Z s D, where the rate of strain is D = (ru + (ru) T )/2. The solvent viscosity Z s is assumed to be constant. In this work the viscoelastic polymer stress is modeled through the constitutive FENE-P model, which is based on the finitely extensible non-linear elastic dumbbell for polymeric materials, as explained in detail by Bird et al. 31 The constitutive equation is derived from kinetic theory, where a polymer chain is represented as a dumbbell consisting of two beads connected by an entropic spring.
Other rheological models, such as the Maxwell model and Oldroyd-B model, take the elastic force between the beads to be proportional to the separation between the beads. The main disadvantage of such models is that the dumbbell can be stretched indefinitely. This leads to a divergent behavior and associated numerical instabilities in strong extensional flows. These problems are prevented by using a finitely extensible spring. The basic form of the FENE-P constitutive equation is: In eqn (3) the operator r above a second-rank tensor represents the upper-convected time derivative, defined as In eqn (3) the constant l is the dominant (longest) relaxation time of the polymer, Z p is the zero-shear rate polymer viscosity, tr(s) denotes the trace of the stress tensor, and L characterizes the maximum polymer extensibility. This parameter equals the maximum length of a FENE dumbbell normalized by its equilibrium length. When L 2 -N the Oldroyd-B model is recovered. The total zero shear rate viscosity of the polymer solution is given as Z = Z s + Z p . The viscosity ratio, which for a real system depends on polymer concentration, is defined as b = Z s /Z. We simulate an unsteady viscoelastic flow through a static array of randomly arranged monodisperse and bidisperse cylinders, constituting a model porous medium, using computational fluid dynamics (CFD). The primitive variables used in the formulation of the CFD model are velocity, pressure and polymer stress. The complete mass and momentum conservation equations are considered in the simulations and discretized in space and time. A coupled finite volume-immersed boundary methodology (FVM-IB) with a Cartesian staggered grid is applied. 1 In the FVM, the computational domain is divided into small control volumes DV and the primitive variables are solved in the control volumes in an integral form over a time interval Dt.
The location of all the primitive variables in a 3D cell is indicated in Fig. 1. The Cartesian velocity components u, v, w are located at the cell faces while pressure p and all components of the stress s are located at the center of the cell. We note that, in this work, we focus on 2D systems where the 3rd dimension has been disabled.
We apply the discrete elastic viscous stress splitting scheme (DEVSS), originally proposed by Guénette and Fortin, 32 to introduce the viscoelastic stress terms in the Navier-Stokes equation because it stabilizes the momentum equation, which is especially important at larger polymer stresses (small b). A uniform grid spacing is used in all directions. The temporal discretization for the momentum eqn (2) is shown as follows, This journal is © The Royal Society of Chemistry 2017 Here Z p r 2 u n+1 and E n p = Z p r 2 u n are the extra variables we introduce to obtain numerical stability, where n indicates the time index. C is the net convective momentum flux given by: Here a first order upwind scheme is used for the implicit evaluation of the convection term (called C f ). In the calculation of the convective term we have implemented a deferred correction method. The deferred correction contribution that is used to achieve second order spatial accuracy, while maintaining stability, is (C n m À C n f ) and is treated explicitly. In this expression C m indicates the convective term evaluated by the total variation diminishing min-mod scheme. A second order central difference scheme is used for the discretization of the diffusive terms.
In eqn (5) the viscoelastic stress part s is calculated by solving eqn (3). The viscoelastic stress tensors are all located in the center of a fluid cell, and interpolated appropriately during the velocity updates. The convective part of eqn (3) is solved by using the higher order upwind scheme.
Eqn (5) is solved by a fractional step method, where the tentative velocity field u** in the first step is computed from: In eqn (7) a set of linear equations is solved. Here it is important to note that the enforcement of a no-slip boundary condition at the surface of the immersed objects is handled at the level of the discretized momentum equations by extrapolating the velocity field along each Cartesian direction towards the body surface using a second order polynomial. 1,33 The main advantage of using the immersed boundary method is that it requires no conformal meshing near the fluid-solid interface whereas the method is computationally robust and cheap. The velocity at the new time step n + 1 is related to the tentative velocity is as follows: where dp = p n+1 À p n is the pressure correction. As u n+1 should satisfy the equation of continuity, the pressure Poisson equation is calculated as: We use a robust and efficient block -incomplete Cholesky conjugate gradient (B-ICCG) algorithm 34,35 to solve the resulting sparse matrix equation for each velocity component in a parallel computational environment.
As the viscoelastic stress tensor components are coupled amongst themselves and with the momentum equation, the velocity at the new time level u n+1 is used to calculate the stress value accordingly. More details and validation of the method can be found in our methodology paper. 1

Problem description
We investigate the flow of a viscoelastic fluid through a static array of randomly arranged cylindrical particles in a 2D periodic domain as shown in Fig. 2. The domain size is set by the solids volume fraction f (or porosity 1 À f), the diameter of each particle d p and number of particles N. To generate the random packing for f r 0.45, a standard hard disk Monte-Carlo (MC) method 36 is employed. The particles are placed initially in an ordered hexagonal configuration in a domain with periodic boundary conditions in all directions. Then each particle is displaced randomly such that no overlap between particles occurs. However, such a MC method does not provide sufficiently random configurations in highly dense packings. 37 Thus, to generate random configurations at f 4 0.45, an event driven method, combined with a particle swelling procedure, is applied. 38 This ensures the particles are randomly distributed. The same approach was followed by Tang et al., for Newtonian fluid simulations for a range of low to intermediate Reynolds numbers. 39 A detailed analysis of different packing generation algorithms and drag correlations for Newtonian fluid flows through such random monodispersed porous media has also been performed. 39,40 In all simulations, the flow is driven by a constant body force exerted on the fluid in the x-direction, while maintaining periodic boundary conditions in x and y directions. Two different pore  geometries are investigated by random arrangement of monodisperse (MD) and bidisperse (BD) sets of cylinders. The diameter ratio of the two cylinders in the bidisperse case is equal to 1.4. 41 The solid fractions f investigated for both bidisperse and monodisperse simulations are 0.3 and 0.6, corresponding to porosities of 0.7 and 0.4, respectively. The viscoelastic fluid is modelled using the FENE-P model, using a constant extensional parameter (L 2 ) of 1000. The viscosity ratio b is kept at 0.33. The objective is to study the interaction between the viscoelastic fluid and solid for different flow configurations and to understand the effect of polydispersity, thus we keep the values of L 2 and b constant. For reference, we also simulate a Newtonian fluid with the same zero-shear viscosity as the polymer solution. In all simulations, we keep the Reynolds number low, below a value of 0.01, ensuring we are always in the creeping flow regime. We perform simulations for Deborah numbers ranging from 0 to 2, where the Deborah numbed is defined as De = lU/R c , based on the polymer relaxation time l, mean flow velocity U, and cylinder radius R c . In the bidisperse system R c is an equivalent radius based on a number average of the areas of the disks, defined as R c ¼ It is important to choose the domain size (or number of particles) sufficiently large to prevent interaction of the flow with its own periodic images. As a consequence, the residence time of the flow in the central periodic domain should be longer than the polymer relaxation time: the length of the domain, in units of R c , should be larger than the maximum Deborah number De. We ensured this was the case in all our simulations. In our previous work, we have investigated finite domain size effects in more detail. 22 To understand the interaction between rheology and porous medium, we will use the concept of flow topology parameter (Q). The flow configuration through the random packings, i.e. the amount of rotational, shear and extensional flow, will depend on the level of viscoelasticity. To characterize the flow configurations, we introduce a flow topology parameter Q which is the second invariant of the normalized velocity gradient. This parameter is defined as where Þ are invariants of the rate of strain tensor D, and the rate of rotation tensor Þ . Values of Q = À1, Q = 0, and Q = 1 correspond to pure rotational flow, pure shear flow and pure elongational flow, respectively. We will further obtain insight in the role of viscoelastic normal stress in such flow configurations by computing the normal stress difference in the frame of reference of local flow and correlating it with the flow topology parameter.

Apparent relative viscosity
The pore configurations, generated by randomly packing monoand bidisperse sets of cylinders, are significantly different for two different solid fractions. Fig. 3 shows viscoelastic flow profiles through a random array of bidisperse cylinders at De = 1.0 for solid fractions (a) f = 0.6 and (b) f = 0.3. As observed from the normalized velocity vectors, the flow pattern is complex in these porous media. For solid volume fraction 0.3, the pore structure is relatively open. Thus there exist regions with less local packing of particles. The fluid prefers to flow through paths connecting these regions of less resistance. However, for solids volume fraction 0.6, the pore structure triggers more tortuous flow paths. Such flow paths enforce the polymers to undergo repetitive contraction and expansion.
To quantify the viscoelastic effects for flow through porous media, we express the results in terms of the viscosity that appears in a generalized Darcy law. The volume-averaged fluid velocity hui in porous media is controlled by the pressure drop across the sample under consideration. According to Darcy's law, for a Newtonian fluid the average pressure gradient (Àdp/dx) and the average fluid velocity across the porous medium are related by: Here Z is the fluid viscosity and k is the permeability (units: m 2 ), which depends on the porosity, pore size distribution and tortuosity of the porous medium. Eqn (11) shows a method to measure the permeability k by flowing a Newtonian fluid of known viscosity through the porous medium. For a viscoelastic fluid, the viscosity is not a constant and generally depends on the flow conditions. However, by assuming that the permeability k is a constant, specific for the porous medium, we can still define an apparent viscosity by using a generalized Darcy law.
Dividing the apparent viscosity by its low flow rate limit  gives us insight in the effective flow-induced thinning or thickening of the fluid in the porous medium. In detail, the apparent relative viscosity Z app of a viscoelastic fluid flowing with a volumetric flow rate q and pressure drop DP through a porous medium is given by: The subscript VE indicates the use of the viscoelastic fluid at a specific flow rate or pressure drop, while the subscript N indicates its Newtonian counterpart in the low flow rate or low pressure drop limit. Fig. 4 shows the change of apparent relative viscosity with an increase in Deborah number (based on cylinder radius) for flow through configurations with different solids volume fractions. For both bidisperse and monodisperse packings, the apparent relative viscosity initially shows a weak shear thinning behavior followed by a flow induced thickening or increased flow resistance after a critical De number. Note that the flow configuration with higher solid volume fraction shows a stronger flow induced resistance (at the same De) and also the onset of increased drag occurs at a lower De number. As the flow configuration with solid fraction of 0.6 has a much stronger fluid solid interaction, the    Up to this point we used the radius of the cylinders as the relevant length scale. However, Fig. 4 suggests that this may not be the appropriate length scale to represent such flow configurations. We next try the square root of the permeability ffiffiffi k p , obtained from Newtonian flow simulations, as the characteristic length scale. This altered Deborah number is defined Fig. 5 shows the apparent relative viscosity versus De k for different solid volume fractions. We find that with this redefined De number, all apparent relative viscosity measurements superimpose to a universal master curve, with onset of flow thickening at De k approximately 1. This is remarkable considering the fact that all geometries have different pore structures, and moreover one system is monodisperse while another is bidisperse. Such a master curve for both mono and bidisperse sets of cylinders, with a high polymer extensibility (high L 2 value) hints that such a master curve might also be obtained for polydisperse or more heterogeneous pore structures.
We have performed simulations of other random configurations at the same porosity. Very similar apparent relative viscosity values (less than 1% difference) were obtained when plotted against De k . This shows that the results can be generalized to any porous medium consisting of randomly placed cylinders. It is important to note that a sufficient amount a randomness in cylinder placement, leading to appearance of preferred flow paths  (as in Fig. 3), is essential for the observed strong increase in apparent relative viscosity. In regular arrays, such as a square array, the increase in apparent relative viscosity is weaker, as found in our simulations (not shown) and in the work of Liu et al. 42 However, even for such regular arrays, we noticed that the onset of this flow thickening starts also at De k approximately 1 (not shown).

Velocity and stress profiles
Next we investigate the profiles (spatial distribution) of velocity and viscoelastic stresses that develop during the flow through the porous media. Although we have investigated two different porosities and both monodisperse and bidisperse sets of cylinders, here we illustrate the profiles for the case of bidisperse cylinders at a solid fraction of f = 0.3. Later we will analyze in detail the interplay between the flow topology and fluid rheology for all cases. Fig. 6 shows snapshots of velocity contours colored by the normalized x-velocity, for different De numbers. With increasing De number, we observe that the flow structure follows more narrowed flow paths with higher velocities, accompanied by more regions with almost stagnant fluid, as shown in Fig. 6(c). Such preferential flow paths emerge due to differences in local flow resistance through various parts of the porous medium, leading to asymmetric flow profiles. We note that in regular periodic arrays of cylinders such asymmetric flow profiles do not occur at the De numbers studied here, 1 although at higher De numbers the symmetry can be broken even in regular arrays. 17 The non-dimensional viscoelastic normal stress component along the average flow direction,t xx /(ZU/R c ), is shown in Fig. 7 for

Flow topology
We next investigate the flow topology. The main purpose is to show how the shear, extensional and rotational parts of the flow are distributed throughout the interstitial pore structure. A flow topology parameter of Q = À1, Q = 0, and Q = +1 corresponds to pure rotational, shear and elongational flow, respectively. Fig. 8 shows the flow topology parameter distribution for a bidisperse random porous medium with solid fraction 0.3. We observe that the flow becomes more shear dominated at higher De. Surprisingly, the presence of extensional flow regions seems to decrease with increasing De.
To show the effect of viscoelasticity on flow topology in a more quantitative way, Fig. 9 shows histograms of flow topology parameter for different De numbers, for each porosity. The common trend, observed for all De numbers, is that the flow structures are more shear dominated than extensional flow dominated. The figure also clearly shows that the prevalence of extensional flow zones definitely does not increase with increasing De number, but rather even decreases (especially clear for the higher solid fraction in Fig. 9(c) and (d)).

Viscoelastic normal stress
The higher prevalence of shear flow zones than extensional flow zones may suggest that shear stresses may be more important in these porous media than extensional stresses. In particular, we suspect that strong normal stress differences in shear dominated regions are responsible for the observed increase in flow resistance at larger De. To investigate this more deeply, we next compute the local first normal stress difference (N 1 ) (normalized by (ZU/R c )) and correlate it with the local flow topology parameter for all locations inside the flow domain.
Here the local first normal stress difference is measured in a local frame of reference, defined by the local flow direction, u/u, and its tangent in the local shear gradient direction. Fig. 10 shows two-dimensional histograms of the mutual occurrence of local first normal stress differences N 1 /(ZU/R c ) and local flow topology parameters Q for a fixed De number. The scale indicated by the color bar is the logarithm of the count. We only show the histograms for the bidisperse flow configurations, but note that the histograms for the monodisperse systems are similar. In case of f = 0.3, at low De number the first normal stresses are distributed quite evenly across both shear (Q = 0) and extensional (Q = 1) dominated flow regimes (Fig. 10(a)). The concentration is actually larger around Q = 1. However, at an increased De of 1.0, we clearly see a concentration of high normal stress difference values around the shear dominated regime around Q = 0. In case of a higher solid fraction of 0.6, at lower De number we can see that the concentration of high normal stress differences is already larger in the shear dominated regime. At higher De number, the maximum normal stress differences are again clearly concentrated around Q = 0. A similar analysis has been performed using the difference between the maximum and minimum eigenvalue of the stress tensor (in the 2D plane). Qualitatively these results are in agreement with Fig. 10 (not shown). These observations strongly suggest that the apparent flow thickening is not primarily due to extensional effects, but rather caused by viscoelastic shear stresses.

Conclusion
We have performed direct numerical simulations using a coupled finite volume immersed boundary methodology to model viscoelastic fluid flow in a 2-dimensional porous medium. The model porous medium is constructed by randomly placing monodisperse and bidisperse sets of cylinders in a periodic domain. Although a shear thinning FENE-P model is employed to model the viscoelastic fluid, we observe an increased flow thickening after a critical De number, irrespective of the pore configuration. The increased flow resistance superimposes to a universal master curve when the characteristic length scale in the definition of Deborah number is selected as the permeability ffiffiffi k p and not the cylinder radius. We have shown a similar master curve in our previous work, 45 for viscoelastic flow through an array of randomly arranged equal-sized spheres representing a three dimensional disordered porous medium. Because the current simulations deal with both mono and bidisperse sets of cylinders, this strengthens our belief that such master curves can also exist for more polydisperse and complex real rock structures, if the length scale is chosen correctly. A detailed analysis of flow topology in the disordered porous media reveals that the flow is more shear dominated than extensional or rotational, especially at higher De number. The analysis of viscoelastic normal stress difference further shows that the largest normal stresses are present in the shear dominated flow regions, rather than in extensional regions. This strongly suggests that viscoelastic normal shear stresses are related to the observed increased flow resistance and not primarily extensional effects. We hope this work has shed new light on the highly debated mechanism governing the increased flow resistance. In our future work, we will study the flow of viscoelastic fluids through realistic rock structures with more complex pore configurations.

Conflicts of interest
There are no conflicts to declare.