Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

S.
De
^{a},
J. A. M.
Kuipers
^{a},
E. A. J. F.
Peters
^{a} and
J. T.
Padding
*^{b}
^{a}Department of Chemical Engineering and Chemistry, Eindhoven University of Technology, The Netherlands
^{b}Process & Energy department, Delft University of Technology, The Netherlands. E-mail: j.t.padding@tudelft.nl

Received
8th September 2017
, Accepted 24th November 2017

First published on 27th November 2017

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.

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–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.

∇·u = 0 | (1) |

(2) |

(3) |

(4) |

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 Δ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 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 β). A uniform grid spacing is used in all directions. The temporal discretization for the momentum eqn (2) is shown as follows,

ρu^{n+1} = ρu^{n} + Δt{−∇p^{n+1} − [C^{n+1}_{f} + (C^{n}_{m} − C^{n}_{f})] + [(η_{s} + η_{p})∇^{2}u^{n+1} + ∇·τ^{n}] + ρg − E^{n}_{p}} | (5) |

C = ρ(∇·uu) | (6) |

In eqn (5) the viscoelastic stress part τ 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:

ρu** = ρu^{n} + Δt{−∇p^{n} − [C_{f}** + (C^{n}_{m} − C^{n}_{f})] + [(η_{s} + η_{p})∇^{2}u** + ∇·τ^{n}] + ρg − E^{n}_{p}} | (7) |

The velocity at the new time step n + 1 is related to the tentative velocity is as follows:

(8) |

(9) |

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}

Fig. 2 Particle configuration at solid fraction ϕ = 0.6 of a random array of monodisperse and bidisperse cylinders (blue represents the interstitial fluid space, the arrow shows the flow direction). |

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 ϕ 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 β 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 β 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 = λU/R_{c}, based on the polymer relaxation time λ, 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 , where N is the total number of cylinders. However here we already note that using the Deborah number De_{k}, based on the length scale equal to the square root of permeability , is probably more important to understand these complex flows. Thus, we will also use De_{k} to analyze different simulation outcomes for different porous configurations.

We have performed simulations for three different mesh sizes: Δ = R_{c}/30, Δ = R_{c}/40 and Δ = R_{c}/80. The results for Δ = R_{c}/40 and Δ = R_{c}/80 were virtually indistinguishable (less than 2% difference in the averaged velocity and stress values), even for De > 1. Therefore, all results in this paper are based on the mesh size Δ = R_{c}/40. We need to keep the CFL number lower than 0.01 in all our simulations, leading to considerable computational costs. At De < 1 a larger time step can be utilized but at De ≥ 1, a smaller time step is required for smooth convergence.

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

(10) |

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.

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 〈u〉 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:

(11) |

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:

(12) |

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 extent of flow resistance is also higher. Even though the bulk fluid is essentially shear thinning in nature, such an increase in flow resistance is consistent with literature findings, especially on experiments in packed beds^{5} and very recent simulations of Liu et al.^{42} Recently, Minale^{43,44} used a generalized Brinkman equation to simulate second order fluid flow through porous media.

Fig. 4 Apparent relative viscosity versus De number for different solid fraction. Here De is based on the (averaged) radius R_{c} of the cylinders. MD = mono-disperse, BD = bi-disperse. |

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 , obtained from Newtonian flow simulations, as the characteristic length scale. This altered Deborah number is defined as . 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).

Fig. 5 Apparent relative viscosity versus De_{k}, using as the characteristic length scale, for different solid fractions. MD = mono-disperse, BD = bi-disperse. |

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,τ_{xx}/(ηU/R_{c}), is shown in Fig. 7 for different De numbers. Such viscoelastic stresses are absent in a Newtonian fluid. Also in our viscoelastic fluids, at a very low De number of 0.01 the dimensionless viscoelastic stresses are very small. With increasing De number the dimensionless viscoelastic stress increases rapidly, where it should be noted that the non-dimensionalization already takes care of the linear increase in flow velocity with increasing De. At a De number of 1.0, we clearly observe the presence of high normalized stresses along the entire flow domain.

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.

Fig. 8 Flow topology parameter Q distribution at (a) De 0.001 and (b) De 1.0, at solid fraction of 0.3. Blue corresponds to rotational flow, green to shear flow, and red to extensional flow. |

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)).

Fig. 9 Flow topology parameter histograms for solid fractions (a) BD ϕ = 0.3, (b) MD ϕ = 0.3, (c) BD ϕ = 0.6, and (d) MD ϕ = 0.6, for different De numbers. |

Fig. 10 shows two-dimensional histograms of the mutual occurrence of local first normal stress differences N_{1}/(ηU/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 ϕ = 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.

- S. De, S. Das, J. A. M. Kuipers, E. A. J. F. Peters and J. T. Padding, J. Non-Newtonian Fluid Mech., 2016, 232, 67–76 CrossRef CAS.
- F. A. L. Dullien, Porous Media-Fluid Transport and Pore Structure, Academic, New York, 1979 Search PubMed.
- L. W. Lake, R. T. Johns, W. R. Rossen and G. A. Pope, Fundamentals of enhanced oil recovery, Society of Petroleum Engineers, 1986 Search PubMed.
- F. A. L. Dullien, Chem. Eng. J., 1975, 10, 1–34 CrossRef CAS.
- R. P. Chhabra, J. Comiti and I. Machao, Chem. Eng. Sci., 2001, 56, 1–27 CrossRef CAS.
- J. G. Savins, Ind. Eng. Chem., 1969, 61, 18–47 CrossRef CAS.
- F. Zami-Pierre, R. De Loubens, M. Quintard and Y. Davit, Phys. Rev. Lett., 2016, 117, 1–5 CrossRef PubMed.
- D. F. James, N. Phan-Thien, M. M. K. Khan, A. N. Beris and S. Pilitsis, J. Non-Newtonian Fluid Mech., 2006, 35, 405–412 CrossRef.
- C. Chmielewski, C. A. Petty and K. Jayaraman, J. Non-Newtonian Fluid Mech., 1990, 35, 309–325 CrossRef CAS.
- R. K. Gupta and T. Sridhar, Rheol. Acta, 1985, 24, 148–151 CrossRef CAS.
- C. Chmielewski and K. Jayaraman, J. Non-Newtonian Fluid Mech., 1993, 48, 285–301 CrossRef CAS.
- K. Arora, R. Sureshkumar and B. Khomami, J. Non-Newtonian Fluid Mech., 2002, 108, 209–226 CrossRef CAS.
- E. S. Shaqfeh, Annu. Rev. Fluid Mech., 1996, 28, 129–185 CrossRef.
- P. Pakdel and G. McKinley, Phys. Rev. Lett., 1996, 77, 2459–2462 CrossRef CAS PubMed.
- A. Groisman and V. Steinberg, Nature, 2000, 405, 53–55 CrossRef CAS PubMed.
- A. Groisman and V. Steinberg, New J. Phys., 2004, 6, 29 CrossRef.
- S. De, J. van der Schaaf, N. G. Deen, J. A. M. Kuipers, E. A. J. F. Peters and J. T. Padding, Phys. Fluids, 2017, 29, 113102 CrossRef.
- L. Pan, A. Morozov, C. Wagner and P. E. Arratia, Phys. Rev. Lett., 2013, 110, 1–5 CrossRef PubMed.
- A. Clarke, A. M. Howe, J. Mitchell, J. Staniland, L. Hawkes and K. Leeper, Soft Matter, 2015, 11, 3536–3541 RSC.
- J. P. Rothstein and G. H. McKinley, J. Non-Newtonian Fluid Mech., 2001, 98, 33 CrossRef CAS.
- S. Louis, J. Fluid Mech., 2009, 631, 231 CrossRef.
- S. De, J. A. M. Kuipers, E. A. J. F. Peters and J. T. Padding, Phys. Rev. Fluids, 2017, 2, 53303 CrossRef.
- D. F. James, R. Yip and I. G. Currie, J. Rheol., 2012, 56, 1249 CrossRef CAS.
- J. P. Singh, S. Padhy, E. S. G. Shaqfeh and D. L. Koch, J. Fluid Mech., 2012, i, 1–37 Search PubMed.
- S. Shahsavari and G. H. McKinley, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 1–27 Search PubMed.
- M. D. Smith, Y. L. Joo, R. C. Armstrong and R. A. Brown, J. Non-Newtonian Fluid Mech., 2002, 109, 13–50 CrossRef.
- A. W. Liu, D. E. Bornside, R. C. Armstrong and R. A. Brown, J. Non-Newtonian Fluid Mech., 1998, 77, 153–190 CrossRef CAS.
- A. Souvaliotis and A. N. Beris, Comput. Methods Appl. Mech. Eng., 1996, 129, 9–28 CrossRef.
- K. K. Talwar and B. Khomami, J. Non-Newtonian Fluid Mech., 1995, 57, 177–202 CrossRef CAS.
- C. C. Hua and J. D. Schieber, J. Rheol., 1998, 42, 477 CrossRef CAS.
- R. B. Bird, R. C. Armstrong and O. Hassager, Dynamics of polymeric liquids, Fluid mechanics, 2nd edn, 1987, vol. 1 Search PubMed.
- R. Guénette and M. Fortin, J. Non-Newtonian Fluid Mech., 1995, 60, 27–52 CrossRef.
- N. G. Deen, S. H. L. Kriebitzsch, M. A. van der Hoef and J. A. M. Kuipers, Chem. Eng. Sci., 2012, 81, 329–344 CrossRef CAS.
- H. A. van der Vorst, Iterative methods for large linear systems, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, 2003 Search PubMed.
- D. S. Kershaw, J. Comput. Phys., 1978, 26, 43–65 CrossRef.
- D. Frenkel and B. Smit, Understanding Molecular Simulation From Algorithms to Applications, 2nd edn, 2001 Search PubMed.
- E. G. Noya, C. Vega and E. De, Miguel, J. Chem. Phys., 2008, 128, 1–7 CrossRef PubMed.
- N. Kumar, O. I. Imole, V. Magnanimo and S. Luding, Adv. Mater. Res., 2012, 508, 160–165 CrossRef.
- Y. Tang, S. H. L. Kriebitzsch, E. A. J. F. Peters, M. A. van der Hoef and J. A. M. Kuipers, Int. J. Multiphase Flow, 2014, 62, 73–86 CrossRef CAS.
- Y. Tang, Direct numerical simulations of hydrodynamics in dense gas-solid flows of Hydrodynamics in Dense Gas-Solid Flows, PhD thesis, TU Eindhoven, 2015 Search PubMed.
- C. S. O. Hern and M. D. Shattuck, Nat. Publ. Gr., 2013, 12, 287–288 Search PubMed.
- H. L. Liu, J. Wang and W. R. Hwang, J. Non-Newtonian Fluid Mech., 2017, 246, 21–30 CrossRef CAS.
- M. Minale, Phys. Fluids, 2016, 28, 023102 CrossRef.
- M. Minale, Phys. Fluids, 2016, 28, 023103 CrossRef.
- S. De, J. A. M. Kuipers, E. A. J. F. Peters and J. T. Padding, J. Non-Newtonian Fluid Mech., 2017, 248, 50–61 CrossRef CAS.

This journal is © The Royal Society of Chemistry 2017 |