Viscoelastic flow simulations in random porous media

We investigate creeping ﬂow of a viscoelastic ﬂuid through a three dimensional random porous medium using computational ﬂuid dynamics. The simulations are performed using a ﬁnite volume methodology with a staggered grid. The no slip boundary condition on the ﬂuid-solid interface is implemented using a second order ﬁnite volume immersed boundary (FVM-IBM) methodology [1] . The viscoelastic ﬂuid is modeled using a FENE-P type model. The simulations reveal a transition from a laminar regime to a nonstationary regime with increasing viscoelasticity. We ﬁnd an increased ﬂow resistance with increase in Deborah number even though shear rheology is shear thinning nature of the ﬂuid. By choosing a length scale based on the permeability of the porous media, a Deborah number can be deﬁned, such that a universal curve for the ﬂow transition is obtained. A study of the ﬂow topology shows how in such disordered porous media shear, extensional and rotational contributions to the ﬂow evolve with increased viscoelasticity. We correlate the ﬂow topology with the dissipation function distribution across the porous domain, and ﬁnd that most of the mechanical energy is dissipated in shear dominated regimes instead, even at high viscoelasticity. open


Introduction
The flow of complex fluids through porous media is a field of considerable interest due to its wide range of practical applications including enhanced oil recovery, blood flow, polymer processing, catalytic polymerization, bioprocessing, geology and many others [2][3][4] . The flow of Newtonian fluids though porous media is relatively well understood in the framework of Darcy's law [2] . Also, a significant effort has been made to understand flow through porous media of non-Newtonian fluids with a viscosity that depends on the instantaneous local shear-rate (inelastic non-Newtonian fluids, or quasi-Newtonian fluids), as reviewed by Chhabra et al. [5] and Savins [6] . However, flow through disordered porous media of viscoelastic fluids, i.e. non-Newtonian fluids displaying elasticity, is far from being understood [5,7,8] . This is due to the complex interplay between the nonlinear fluid rheology and the porous geometry. Several types of numerical frameworks have been used to model flow of non-Newtonian fluids through porous media, including extensions of Darcy's law [9] , capillary based models [10] , and direct numerical simulations based on computational fluid dynamics. Unfortunately, extensions of Darcy's law and capillary based models are found to be inadequate to * Corresponding author. accurately capture the complete physics of pore scale viscoelastic flow through porous media [11][12][13] .
Many numerical works focus on relatively simple geometries to uncover the essentials of non-Newtonian fluid flow through porous media [14][15][16][17] . Sometimes a full three-dimensional random porous medium is studied, which is already closer to a realistic pore geometry, but such studies are then usually limited to power-law fluids, which are the most commonly applied quasi-Newtonian fluids [11,[18][19][20] . For example, Morais et al. [18] applied direct numerical simulations to investigate the flow of power-law fluids through a disordered porous medium and found that pore geometry and fluid rheology are responsible for an increase in hydraulic conductance at moderate Reynolds numbers. Simulations of fully viscoelastic fluid flows are limited to two dimensional pore geometries [21][22][23][24][25] . It is now commonly agreed that including viscoelasticity is important: both numerically and experimentally, viscoelasticity is found to introduce profound effects and complex phenomena such as enhanced pressure drop and elastic instabilities (sometimes referred to as elastic turbulence) [5,[26][27][28][29][30][31][32][33][34][35] . So, although it is known that viscoelastic fluids behave more complex than inelastic non-Newtonian fluids, the current literature shows a lack of detailed simulations of fully three dimensional flows of viscoelastic fluids through random porous media.
In this paper, we report on a numerical study of the flow of viscoelastic fluids through three dimensional random porous media consisting of packed arrangements of monodispersed spherical particles using a combined finite volume immersed boundary (FVM-IBM) methodology. Four different porosities are studied for a range of low to high Deborah numbers (defined later). We measure in detail the viscoelastic fluid flow structure and stress development in the porous medium. We will show a transition from a symmetric Newtonian flow profile to an asymmetric flow configuration, and will relate it to a strong increase in pressure drop. An analysis of the flow topology will show how shear, extension and rotation dominated flow regimes change with increasing viscoelasticity for different porous structures. Finally, we will show how the distribution of mechanical energy dissipation in the porous medium changes with increasing viscoelasticity and correlate this with the flow topology. This analysis will help us to understand the interplay of pore structure and fluid rheology in three dimensional random porous media.

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, ρ is the fluid density (assumed to be constant) and p is the pressure. τ is the viscoelastic polymer stress tensor. The Newtonian solvent contribution is explicitly added to the stress and defined as 2 η s D , where the rate of strain is D = (∇u + (∇u ) T ) / 2 . The solvent viscosity η 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. [36] . The equation is derived from a kinetic theory, where a polymer chain is represented as a dumbbell consisting of two beads connected by an entropic spring. Other basic 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, leading to 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 Eq. (3) the operator ∇ above a second-rank tensor represents the upper-convected time derivative, defined as In Eq. (3) the constant λ is the dominant (longest) relaxation time of the polymer, η p is the zero-shear rate polymer viscosity, tr ( τ) 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 → ∞ the Oldroyd-B model is recovered. The total zero shear rate viscosity of the polymer solution is given as η = η s + η p .
The viscosity ratio, which for a real system depends on polymer concentration, is defined as β = η s /η.
We simulate an unsteady viscoelastic flow through a static array of randomly arranged monodisperse spheres, constituting a model porous medium, using computational fluid dynamics (CFD). The primitive variables used in the formulation of the model are velocity, pressure and polymer stress. The complete mass and momentum conservation equations are considered and discretized in space and time. A coupled finite volume -immersed boundary methodology [1] (FVM -IBM) with a Cartesian staggered grid is applied. In the FVM, the computational domain is divided into small control volumes V and the primitive variables are solved in the control volumes in an integral form over a time interval t .
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 τ are located at the center of the cell.
We apply the discrete elastic viscous stress splitting scheme (DEVSS), originally proposed by Guénette and Fortin [37] , to introduce the viscoelastic stress terms in the Navier-Stokes equation because it stabilizes the momentum equation, which is especially im- Here η p ∇ 2 u n +1 and E n p = η p ∇ 2 u n are the extra variables we introduce to obtain numerical stability, where n indicates the time index. C represents the net convective momentum flux given by: Here the 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 (CD) scheme is used for the discretization of diffusive terms. Eq. (5) is solved by a fractional step method, where the tentative velocity field in the first step is computed from: In Eq. (7) we need to solve a set of linear equations. 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,38] . 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 δ p = 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 [39,40] to solve the resulting sparse matrix 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. As a steady state criterion, the relative change of velocity and stress components between two subsequent time steps are computed in all the cells. If the magnitude of the relative change is less than 10 −4 the simulation is stopped. This part is also explained in detail in our methodology paper [1] .

Problem description
We employ our method to investigate the flow of viscoelastic fluid through a static array of randomly arranged spherical particles in a 3D periodic domain ( Fig. 2 ). The domain size is set by the solids volume fraction φ, the diameter of each particle d p and number of particles N p . To generate the random packing for φ ≤ 0.45, a standard hard sphere Monte-Carlo (MC) method [41] is used. The particles are placed initially in an ordered face centered cubic (FCC) configuration in a domain with periodic boundary conditions in all directions. Then each particle is moved randomly such that no overlap between particles occurs. However, such a MC method does not provide sufficiently random configurations in highly dense packings [42] . Thus, to generate random configurations at φ > 0.45, an event driven method combined with a particle swelling procedure is applied [43] . 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 [44] . A detailed analysis of different packing generation and drag correlation study for Newtonian fluid flow through such a random monodispersed porous media has also been performed [44,45] .
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 all three directions. Note that this is slightly different from using periodic boundary conditions with a pressure jump condition. Simulations of random arrays are carried out with N p = 108 spheres arranged in different configurations. The  The precise 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 configuration, we introduce a flow topology parameter Q which is the second invariant of the normalized velocity gradient. This parameter is defined as where S 2 = 1 2 ( D : D ) and 2 = 1 2 ( : ) are invariants of the rate of strain tensor D , introduced before, and the rate of rotation tensor = 1 2 ( ∇ u T − ∇u ) . Values of Q = − 1, Q = 0, and Q = 1 correspond to pure rotational flow, pure shear flow and pure elongational flow, respectively. In this paper we will correlate the above flow topology parameter Q with the dissipation function in the flow domain. The dissipation function expresses the work performed by the viscoelastic and Newtonian stress per unit volume (in W/m 3 ), and is defined as: . e t = ( τ + 2 η s D ) : ∇u (11) By correlating the spatial distributions of Q and . e in the porous domains at different De numbers, we will be able to identify the flow configurations which lead to the predominant energy dissipation, and are therefore predominantly responsible for observed pressure drops.
To quantify the dissipation function in a dimensionless manner, we express the total work performed by Newtonian and viscoelastic stress per unit of volume as E t = . e t ηU 2 /R 2 C . We warn the reader that the word "dissipation function" may be a bit misleading because for a viscoelastic fluid not all of the work represented by this term is irreversibly turned into heat, but instead can be stored elastically and released at a later point in time, leading to a local negative value of the dissipation function. We will show later that this indeed is the case. These streamlines provide an idea about the complex flow pattern in these porous media. For solids volume fraction 0.3, the flow is rather homogeneous. However for solids volume fraction 0.5, the pore structure triggers more tortuous flow paths and more preferential flows.

Apparent relative viscosity
To quantify the viscoelastic effects we express the results in terms of the viscosity that appears in a generalized Darcy law for flow through porous media. The volume-averaged fluid velocity u in porous media is controlled by the pressure drop across the sample. According to Darcy's law, for a Newtonian fluid the relation between the average pressure gradient ( − dp dx ) and the average fluid velocity across the porous medium is: Here k is the permeability (units: m 2 ), which is related to the solids volume fraction (or porosity), pore size distribution and tortuosity of the porous medium, whereas η is the viscosity of the fluid. Eq. (12) presents an operational way of measuring 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 but generally depends on the flow conditions. However, if we assume that k is constant for a specific 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 η app of a viscoelastic fluid flowing with a volumetric flow rate q and pressure drop P through a porous medium is given by: The subscript VE indicates 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 depicts how the apparent relative viscosity changes with an increase in viscoelasticity for flow through flow configurations with different solids volume fractions. With increasing De number, where De is based on the sphere radius as the characteristic length scale, we initially observe a (relatively weak) flow-induced thinning. Then beyond a certain flow rate we observe a strong flow-induced thickening, which means a sharp increase in flow resistance. With increasing solids volume fraction (decreasing poros-ity), the onset of this increased flow resistance shifts to a lower De number. This shows that the increased fluid-solid interaction facilitates the onset of such a flow resistance. Experimental evidence of this increase in apparent relative viscosity was previously reported in literature [5] , especially for packed bed systems.
The pore porosity and pore geometry are very important for the increase in apparent relative viscosity, but this is not reflected in the De number based on the radius of the spheres. Therefore, we next try to use the square root of the permeability, √ k obtained from Newtonian flow simulations, as the characteristic length scale. This altered Deborah number is defined as D e k = λU √ k . Fig. 5 shows the apparent relative viscosity versus the altered De k for different solids volume fractions. We find a collapse of all data sets of Fig. 4 to a single curve for the entire range of De k numbers. This is remarkable considering the fact that, despite the different arrangement of pore structures for the different porosities, the resulting increase in flow resistance follows the same universal thickening behavior. However we should keep in mind that these results are strictly only valid for a FENE-P type of fluid with L 2 = 100 flowing through a random array of monodisperse spheres.  The increase in flow resistance for flow of viscoelastic fluid through packed bed are also experimentally shown in the work of Chhabra et al. [5] and by W. Kozicki [46] based on a capillary hybrid model. Recently M. Minale, [47,48] showed a similar kind of scaling relationship for flow of second order fluids through porous media using a generalized Brinkman's equation.

Velocity and stress profiles
Next we investigate the velocity and stress profiles of viscoelastic fluid flow through the three dimensional porous medium, and analyze the interplay between the flow structures and fluid rheology. Although we have investigated different porosities, here we show the profiles for a solid fraction of φ= 0.5 for a range of De numbers. Fig. 6 shows snapshots of velocity contours (across a representative section of flow domain), colored by the normalized xvelocity, for different De numbers after the same time of simulation. The flow structure becomes non uniform with increasing De number. Especially at De of order 1 we see the onset of preferential flow paths (paths with higher velocity) in the flow domain. Fig. 7 illustrates the same effect with streamlines (colored with normalized vorticity), clearly showing the meandering flow paths through the pore space. Such preferential flow paths emerge due to differences in flow resistance through different parts of the porous medium, leading to asymmetric flow structures, as will be discussed in detail later.
The non-dimensional viscoelastic normal stress component along the flow direction τxx ηU/ R c , is shown in Fig. 8 for different De numbers. Such viscoelastic stresses are absent in a Newtonian fluid. We observe that the viscoelastic normal stress increases with increasing De number, and that the largest normal stresses are present near the walls of the sphere at locations which are shear dominated. This will also be analyzed in detail in the subsequent section.
Though the flow is in a non-inertial regime, at higher viscoelasticity the uniformity in the streamlines is found to be less, showing flow asymmetry, compared to its Newtonian counterpart.
To understand the effect of viscoelasticity on flow anisotropy we have analyzed the velocity probability distribution function (PDF) across the entire three dimensional flow domain (each porosity) for all the De numbers based on √ k . Fig. 9 shows the distribution of normalized velocities along the flow direction x. For low Reynolds numbers such as studied here one might expect the PDFs to collapse on each other. This is not the case. At low De number the PDFs of the x velocity component superimpose and are mostly positive. However, at increased De numbers the PDFs also increase for negative velocities. This shows that there is emergence of recirculation zones in the system. Though the driving force for the flow is along the positive x direction, the negative components in the PDFs give a measure for recirculation appearing in the system. Another interesting fact is that the width of the PDFs increases with increasing solids volume fraction, and show an appearance of a slower decaying tail for higher velocities. Fig. 10 shows the distribution of normalized velocities along the transverse flow direction y. In a non-inertial flow regime with random placement of the spheres we expect a symmetric distribution of y-velocities. Fig. 10 , shows that for low De number the PDFs of the transverse velocity components are completely symmetric. However with increasing De number the PDFs become slightly asymmetric, and the broadness of the PDFs also increases. The decay is more exponential in nature. The reader should keep in mind that the vertical scales of the PDF distributions are plotted on a log scale, so the tails represent very small probabilities. Thus, because a finite number of samples are used, the PDFs are not smooth in the tails.
These observations quantitatively validate our findings of the existence of preferential flow paths, observed from the streamlines. A possible mechanism might be that, at increased viscoelasticity, strong elastic effects come into play, leading to asymmetric curved streamlines and possibly causing elastic instabilities, as also shown in the work of Pakdel et al. [30] . To understand these effects further we have performed a detailed analysis of the flow topology and dissipation function, and will be presented in the next section.

Flow topology
This section focuses on the flow topology. As explained in Section 2.2 , the main idea is to investigate how the shear, extensional and rotational parts of the flow are distributed and develop in the three dimensional interstitial space. As explained Q = −1, Q = 0, and Q = + 1 correspond to pure rotational, shear and elongational flows, respectively. Fig. 11 shows the flow topology parameter distribution for a random porous medium with solid fraction 0.5. We observe that the flow becomes more shear dominated at higher De, while, perhaps surprisingly, the presence of extensional flow regions seems to decrease. To better quantify the effect of viscoelasticity on flow topology, in Fig. 12 we have plotted the histograms of flow topology parameter for different De numbers, for each solids volume fraction. The common feature observed from all histograms is that all flow structures are more shear dominated than extensional flow dominated. Although the extensional component ( Q = 1) increases slightly up to De = O(1) , at larger De it sharply decreases again and shear effects ( Q = 0) become more dominant. Note that the PDFs of normalized velocity ( Fig. 10 ) and flow streamlines ( Fig. 7 ) show a transition to flow asymmetry around the same De = O(1) . So the flow topology analysis shows that the increase in flow resistance at larger De, observed in a random porous medium of monodisperse spheres, may be caused by strong normal viscoelastic shear flow stresses, rather than their extensional counterparts. Fig. 13 compares the flow topology histograms at the same De k = 1.0 for four different solids volume fractions. This shows that at this relatively high De number, the overall shear contribution (Q = 0) also increases with increasing solids volume fraction (decreasing porosity), and subsequently the extensional contribution decreases. These results are also consistent with our recent observations of viscoelastic fluid flow in a model porous media [49] .

Dissipation function
Finally, we analyze the spatial distribution of the dissipation function, expressing the work done by the total stress (viscoelastic + Newtonian solvent) per unit of time and per unit of volume, as defined in Section 2.2 . This dissipation function can be both positive and negative, but we note that energy is always dissipated from the Newtonian solvent contribution. As an example, Fig. 14 shows the nondimensionalised spatial distribution of E t for a solid fraction of 0.5 and De number of 1.0, in a representative plane of the random porous media. We have clipped the dissipation function color scale to clearly show regions in the domain where energy is released (negative dissipation function) by the polymer solution. For De = O(1) and higher, energy is dissipated by the solvent, but also stored as elastic energy by the polymers, close to the particle surfaces, and released after a pore throat has ended, further away from the particle surfaces. This is consistent with the physical picture in which polymers in fast contraction flow are extended and therefore store energy in their entropic springs; this energy is subsequently released when the polymers can recoil to their upstretched state when the contraction flow has stopped.
In the previous section, we showed that the fraction of shear dominated regions increases significantly beyond De = O(1) . We now ask whether these shear dominated regions are also responsible for the observed increase in flow resistance. To answer this, in Fig. 15 we show what fraction of energy is dissipated in the flow domain with a particular value of flow topology Q (per unit Q ).
Only for the lowest solid fraction 0.3 (highest porosity of 0.7), we find significant energy dissipation in mixed shear and extensional flow (0 < Q < 1). For more closely packed domains this zone significantly reduces: the width of the histograms reduce and their peaks grow around Q = 0 with increasing solid fraction, while the contribution of extensional flow to the energy dissipation generally decreases with increasing De. This conclusively shows that at increased solid fraction and with increasing De, shear regions are predominantly responsible for the increase in flow resistance in the random porous media studied here.

Conclusion
We have employed a finite volume -immersed boundary methodology to study the flow of viscoelastic fluids through an array of randomly arranged equal-sized spheres representing a three dimensional disordered porous medium, for a range of solid fractions (or porosities). Irrespective of the solid fraction, we found a strong increase in flow resistance after a critical De number is reached. The increase in apparent relative viscosities measured for different solids volume fractions overlap among each other if the Deborah number is chosen with a length scale based on the permeability of the pore space (more precisely, D e k = λU/ √ k , with k the permeability of the medium for a Newtonian fluid. The PDFs of flow velocity suggest that with increasing viscoelasticity the flow profiles become more asymmetric, and increasingly preferential flow paths are found. A detailed study of the flow topology shows that for the porous media investigated in our study, shear flow becomes more important than extensional or rotational flow at higher De number. We have analyzed the distribution of the dissipation function across the flow domain and correlated it to the flow topology. These findings helped us conclude that the observed increase in flow resistance should be attributed to an increase in energy dissipation in shear flow dominated regions. More generally, simulations such as shown here help us to understand the complex interplay between the fluid rheology and pore structure in porous media. In our future work we will study flow through three dimensional realistic porous media which have a larger distribution in pore and throat sizes than studied here.