Elastoviscoplastic flows in porous media

https://doi.org/10.1016/j.jnnfm.2018.04.006Get rights and content

Highlights

  • Shear and elongational flow are equally distributed in the unyielded region.

  • The apparent permeability decreases with the Bingham number.

  • The unyielded region increases with the Weissenberg number.

  • Relation between pressure drop, flow rate and Bingham number.

Abstract

We investigate the elastoviscoplastic flow through porous media by numerical simulations. We solve the Navier–Stokes equations combined with the elastoviscoplastic model proposed by Saramito for the stress tensor evolution [1]. In this model, the material behaves as a viscoelastic solid when unyielded, and as a viscoelastic Oldroyd-B fluid for stresses higher than the yield stress. The porous media is made of a symmetric array of cylinders, and we solve the flow in one periodic cell. We find that the solution is time-dependent even at low Reynolds numbers as we observe oscillations in time of the unyielded region especially at high Bingham numbers. The volume of the unyielded region slightly decreases with the Reynolds number and strongly increases with the Bingham number; up to 70% of the total volume is unyielded for the highest Bingham numbers considered here. The flow is mainly shear dominated in the yielded region, while shear and elongational flow are equally distributed in the unyielded region. We compute the relation between the pressure drop and the flow rate in the porous medium and present an empirical closure as function of the Bingham and Reynolds numbers. The apparent permeability, normalized with the case of Newtonian fluids, is shown to be greater than 1 at low Bingham numbers, corresponding to lower pressure drops due to the flow elasticity, and smaller than 1 for high Bingham numbers, indicating larger dissipation in the flow owing to the presence of the yielded regions. Finally we investigate the effect of the Weissenberg number on the distribution of the unyielded regions and on the pressure gradient.

Introduction

Fluid flows through porous media are relevant for different industrial applications such as oil recovery, polymer extrusion, filtration processes and food processing. They are also present in biological flows that involve mass transfer across organic tissues, like blood vessels, kidney and lungs. Sedimentary rocks and riverbeds are examples of this kind of flows in nature [2]. The difficulty in understanding these flows arises from the complex structure of the porous media as well as from its multiscale nature; often, the non-Newtonian behaviour of the involved fluids further complicates the dynamics.

The objective of this work is to investigate the elastoviscoplastic (EVP) flow through a porous medium and understand the role of yield stress in the relationship between the pressure drop and the flow rate. Hence, the focus of the work will be on non-Newtonian effects in porous media.

Porous media can be defined as a solid matrix with interconnected cavities distributed inside. Usually, these materials are described from a macroscopic point of view, yet the relevant parameters strongly depend on the microscopic structure. The characteristic properties of a porous medium are the porosity and the permeability. The porosity, ε, measures the void space inside the material and is defined as the ratio between the volume of the voids and the total volume of the medium, thus, varying between 0 and 1. The permeability, K, measures the ability of the flow to pass through the porous medium and has the dimension of a length squared. If the medium is impermeable K=0 whereas, if the medium offers no resistance to a flow, the permeability is infinite. Typical values of porosity for porous media can range from 0.32 for cylindrical packings, to 0.8 for foam metals [3], [4]. Most synthetic and natural porous media exhibit inhomogeneity due to the randomness of their structure. However, in literature is a common practice to consider as a first approximation the medium as homogeneous and composed of an array of cylinders and consider only a periodic cell [5], [6], [7], thus we adopted the same discretisation in this study.

From an engineering point of view, it is of primary importance to determine the pressure gradient required to obtain a desired flow rate through the porous medium. The pioneer in the field was Darcy who, in 1856, derived an empirical one-dimensional relation between flow rate and pressure drop, i.e., the well-known Darcy law [8]. The law was then extended into three dimensions, see for example [9], and also derived theoretically by Whitaker [10]. Note that this analytical relation between the pressure drop and the flow rate across the porous media is valid only for Newtonian fluids. An extension of Darcy law to non-Newtonian flows has been derived using homogenisation theory, see [11]. In particular, Hornung [11] states that the system of equations derived with the homogenisation method reduces to a non-linear Darcy-like law only if the flow is directional and that “We can not always expect to have laws as simple as Darcy’s, describing complex natural phenomena”. It is worth noticing that the law derived in Ref. [11] is based on the assumption of a steady-state flow and in the absence of elasticity effects, i.e. generalized Newtonian inelastic fluids. However, many real flows in porous media exhibit both elastic and plastic properties. Additionally, other effects that produce non-Darcy behaviours should be considered, inertia being one example of those; see for example Refs. [12], [13].

The fluids flowing inside porous media exhibit often a non-Newtonian behaviour, characterised by strain- and time-dependent viscosity, yield-stress and/or stress relaxation [14]. Many models (constitutive equations) have been proposed in literature to describe these materials, from simple power laws for inelastic shear-dependent viscosity fluids to more complex visco-elastic models like the Oldroyd-B closure. The reader is refereed to Ref. [15] for a detailed review of non-Newtonian flows in porous media.

Viscoelastic fluids exhibit both viscous-fluid and elastic-solid behaviour. Polymer solutions are often modelled as viscoelastic materials. Their behaviour can be described combining the Newton law for viscous fluid with Hooke law for elastic solid, as done originally by Maxwell [16]. In viscoelastic fluids, the stress depends on the history of the rate of strain, with a fading memory. Hence, for slow deformation viscoelastic materials behave as viscous fluids whereas for fast deformation they behave as elastic solids. For this reason, the important parameter characterising the flow of viscoelastic fluids is the ratio between the material time scale and the time scale of the flow, indicated by the Deborah number (De) [17].

Flow of viscoelastic materials through porous media has been the object of numerous research in the past owing to the large range of applications involved. Experimental studies have been performed to investigate the pressure drop changes due to elasticity [18], [19], [20] or the instability onset for high Deborah numbers [5], [21], [22], [23]. Numerical simulations have also proved to be a valid tool to investigate viscoelastic flows in porous media, uncovering the details of the fluid deformation and stresses. Alcocer and Singh [24] studied the viscoelastic flow through an array of cylinders using the finite element method and a FENE (finitely extensible nonlinear elastic) model for the stresses. These authors have shown the influence of the cylinder distribution and Deborah number on the permeability of the medium. Richter et al. [25] studied the effects of the Reynolds number on the instability of a FENE-P viscoelastic fluid flowing around a cylinder and reported an increase in the drag for higher extensibility of the polymer as well as the suppression of the three-dimensional Newtonian B-mode instability for Re ∼ 300. A similar study, but for the flow around a sphere is presented in Ref. [26]. The effect of high-Weissenberg number (defined as the ratio between elastic force and viscous force) has been investigated for the Oldroyd-B and Giesekus model in Ref. [27], [28]. It is shown that by mean of the log conformation representation it is possible to obtain average convergence of the solution beyond the limiting Weissenberg. More recently, Grilli et al. [29] have simulated the elastic turbulent flow of an Oldroyd-B fluid through an array of cylinders. De et al. [7] have anaylzed the effect of the arrangement of the cylinders on the flow at low Reynolds numbers. They have found differences between symmetric and asymmetric configuration in the flow characteristic and effective viscosity. The effects of elasticity on a flow through a porous medium composed of a random distribution of spheres has been studied in [30], [31] while the effect of the cylinder distribution on the elastic instability has been investigated in [32]. Similarly, Shahsavari and McKinley [33] studied the flow of a Carreau fluid through a periodic array of cylinders and proposed a scaling for the mobility functions.

As already stated before, viscoelasticity is not the only property exhibited by non-Newtonian fluids. Another important characteristic is the yield stress or the ability to sustain shear stress. A yield stress fluid, or visco-plastic fluid, behaves as a solid below the yield stress and as a fluid above. A list of several materials exhibiting yield stress is provided in [34]. The first attempt to study the rheology of a viscoplastic material dates back to Schwedoff [35], based on the Maxwell model. Next, Bingham [36] proposed a one-dimensional stress-deformation rate equation that was extended into three dimensions and coupled with the Navier–Stokes equations by Oldroyd [37]. The latter also proposed a new constitutive equation considering a linear Hookean behaviour before yielding. According to this model, the material is no longer rigid before yielding and the stress is expressed as function of the strain and not the strain rate, leading to a cuspid in the constitutive equation.

Experiments of sediments and particles interactions in Carbopol solutions and Laponite suspensions [38], [39], [40] have shown the loss of symmetry corresponding to the elastic effects. Therefore, it is essential to include elastic effects in dealing with conventional yield stress test fluids such as Carbopol, yet only few works have been conducted to introduce elasticity alongside plasticity.

Recently, Saramito [1] proposed a constitutive equation for EVP fluid flows based on thermodynamic principles. This model reproduces a viscoelastic solid for stresses lower then the yield stress whereas the material behaves as a viscoelastic Oldroyd-B fluid for stresses higher then the yield stress. To describe the fluidisation process it uses the von Mises yielding criterion, which has been experimentally confirmed [41], [42]. Other models based on the Papanastasiou regularisation have also been proposed to model yield stress [43], [44], see [45] for a detailed analysis of these models.

The model proposed in [1] was extended by the same author to account for shear-thinning effects [46]. The new model combines the Oldroyd-B viscoelastic model with the Herschel–Bulkley viscoplastic model, a power law with index n > 0. When n=1 the model reduces to the one proposed in [1]. Benito et al. [47] derived a minimal, fully tensorial, rheological constitutive equation for EVP flows to study the behaviour of immortal fluids: it can describe a large deformation in the elastic regime and predict viscoelastic fluid after yielding. To avoid the zero loss modulus at low shear, Dimitriou et al. [48], [49] proposed a new constitutive model based on the Isotropic Kinematic Hardening (IKH) idea. This model can predict thixotropic behaviours and has been proven to correctly describe EVP materials such as waxy crude oils. Armstrong et al. [50] modified the Delaware thixotropic model and proposed a new structure-based model which accounts for shear-induced aggregation. This new model has been shown to correctly describe SAOS (small amplitude oscillatory shear) and weak LAOS (large amplitude oscillatory shear) flows. More recently, Wei et al. [51] proposed a new model to predict transient thixotropic EVP (TEVP) flow. They combined the multi-lambda (ML) model of Wei et al. [52] with the IKH of Dimitriou and McKinley [49] deriving a new constitutive equation, called ML-IKH, with 12 parameters. This new model accounts for a non-linear thixotropic kinetic equation, kinematic hardening and linear viscoelasticity and has been proven to correctly describe TEVP flows in many conditions, such as LAOS and shear reversal.

Many attempts have been carried out in recent years to derive experimentally a relation between the pressure drop and the flow rate in porous media for a yield stress fluid [53], [54], [55], [56]. Balhoff and Thompson [57] noted that it is not clear if a general relation can be found for yield stress fluid flows or it should be derived case-by-case. Cheddadi et al. [58], in particular, have studied the flow of an EVP fluid around an obstacle and shown that viscous, elastic and plastic effects cannot be considered separately before and after yielding but need to be taken into consideration simultaneously. Only few numerical studies have been conducted for yield-stress and, in particular, for EVP flows, see [59] and references therein. Recently, Roustaei et al. [13] have investigated numerically the flow of a yield stress fluid along an uneven fracture. They have shown that the approximation error consequent to the evaluation of the pressure drop using a Darcy-type law strongly depends on the geometry of the problem.

We therefore simulate the EVP flow through a model porous medium, represented by an array of cylinders; here only a symmetric periodic cell is considered. The non-Newtonian flow is simulated by solving the full incompressible Navier–Stokes equations coupled with the model proposed by Saramito [1] for the evolution of the EVP stress tensor. In the next section, we present the governing equations and their numerical discretisation, as well as the validation performed to verify our mathematical formulation and its numerical implementation. We describe the chosen flow configuration in Section 3 and discuss the results in Section 4. In particular, we examine the influence of the Reynolds and Bingham numbers on the flow in the porous media. Finally, we summarise the main findings and conclusions in Section 5.

Section snippets

Formulation

The dynamics of an incompressible EVP flow is fully described by the Navier–Stokes equations with an additional constitutive equation for the non-Newtonian stress tensor. The Navier–Stokes equations expressing the mass conservation and momentum balance read·u=0,ρDuDt=p+·(2μfD)+ρf+·τ.

In the previous set of equations, the symbol D/Dt=/t+u· denotes the material derivative sum of the time derivative and the advection term, u=u(x,t) and p=p(x,t) are the velocity and pressure fields, ρ and μf

Problem description

We consider the incompressible EVP flow through a model porous medium composed by an array of cylinders with porosity ɛ=0.38. Here, we focus on a single periodic cell, as sketched in Fig. 5, following the work by [7].

The numerical domain is a square box of size L=2.25r, r being the radius of the cylinders centered in each corner of the domain and the length-scale of the problem. A periodic boundary condition is enforced in the streamwise xdirection, while the free-slip boundary condition is

Results

First, we provide a qualitative picture of the flow through the porous medium considered hereafter. Fig. 6 shows streamlines and contours of the streamwise velocity component at a Reynolds number of 0.01 for a Newtonian fluid. The flow enters the domain in the narrow neck between two adjacent cylinders, strongly decelerates due to the sudden expansion of the geometry and then accelerates again due to the contraction formed by the following series of cylinders. Close to the top and bottom

Conclusion

We have performed numerical simulations of the elastoviscoplastic flow through porous media modelled as an array of cylinders, and considered a single periodic cell. The flow is described by the Navier–Stokes equations, and the additional evolution equation for the EVP stress tensor following the model proposed by Saramito [1].

We find the flow dynamics to be time-dependent, and have thus presented both instantaneous configurations and time-averaged results. The data show that the volume where

Acknowledgements

This work was supported by the European Research Council Grant no. ERC-2013-CoG-616186 and TRITOS. We also acknowledge financial support by the Swedish Research Council through grants no. VR 2013-5789, no. VR 2014-5001 and no. VR 2017-76478. S.H. acknowledges financial support by NSF (Grant no. CBET-1554044-CAREER), NSF-ERC (Grant no. CBET-1554044 Supplementary CAREER) and ACS PRF (Grant no. 55661-DNI9). The authors acknowledge computer time provided by SNIC (Swedish National Infrastructure for

References (77)

  • S. De et al.

    Viscoelastic flow simulations in random porous media

    J. Non-Newton. Fluid Mech.

    (2017)
  • H.L. Liu et al.

    Flow resistance of viscoelastic flows in fibrous porous media

    J. Non-Newton. Fluid Mech.

    (2017)
  • D. Fraggedakis et al.

    Yielding the yield stress analysis: a thorough comparison of recently proposed elasto-visco-plastic (evp) fluid models

    J. Non-Newton. Fluid Mech.

    (2016)
  • P. Saramito

    A new elastoviscoplastic model based on the herschel-bulkley viscoplastic model

    J. Non-Newton. Fluid Mech.

    (2009)
  • T. Chevalier et al.

    Darcy’S law for yield stress fluid flowing through a porous medium

    J. Non-Newton Fluid Mech.

    (2013)
  • E. Fadlun et al.

    Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations

    J. Comput. Phys.

    (2000)
  • J. Kim et al.

    Application of a fractional-step method to incompressible navier-stokes equations

    J. Comp. Phys.

    (1985)
  • R. Fattal et al.

    Time-dependent simulation of viscoelastic flows at high weissenberg number using the log-conformation representation

    J. Non-Newton. Fluid Mech.

    (2005)
  • D. Izbassarov et al.

    A front-tracking method for computational modeling of viscoelastic two-phase flow systems

    J. Non-Newton. Fluid Mech.

    (2015)
  • S. Hormozi et al.

    Stable core-annular flows of viscoelastic fluids using the visco-plastic lubrication technique

    J. Non-Newton. Fluid Mech.

    (2011)
  • S. Hormozi et al.

    Nonlinear stability of a visco-plastically lubricated viscoelastic fluid flow

    J. Non-Newton. Fluid Mech.

    (2012)
  • S. Hormozi et al.

    Visco-plastic sculpting

    Phys. Fluids

    (2014)
  • F.A. Dullien

    Porous Media: Fluid Transport and Pore Structure

    (2012)
  • I.F. Macdonald et al.

    Flow through porous media-the ergun equation revisited

    Indus. Eng. Chem. Fundam.

    (1979)
  • G.S. Beavers et al.

    Boundary conditions at a naturally permeable wall

    J. Fluid Mech.

    (1967)
  • S. De et al.

    Viscoelastic flow simulations in model porous media

    Phys. Rev. Fluids

    (2017)
  • H. Darcy

    Les Fontaines Publiques de la Ville de dijon

    Dalmont, Paris

    (1856)
  • R.A. Greenkorn

    Flow Phenomena in Porous Media: Fundamentals and Applications in Petroleum, Water, and Food Production

    (1984)
  • S. Whitaker

    Flow in porous media i: A theoretical derivation of darcy’s law

    Transp. Porous Media

    (1986)
  • U. Hornung

    Homogenization and Porous Media

    (2012)
  • S. Kristlanovitch

    Motion of ground water which does not conform to darcy’s law

    Prils. Mat. Mech. (Russian)

    (1940)
  • A. Roustaei et al.

    Non-darcy effects in fracture flows of a yield stress fluid

    J. Fluid Mech.

    (2016)
  • L.R. G.

    Constitutive Equations for Polymer Melts and Solutions

    (1988)
  • R.J. Poole

    The Deborah and Weissenberg numbers

    Rheol. Bull.

    (2012)
  • R.J. Marshall et al.

    Flow of viscoelastic fluids through porous media

    Indus. Eng. Chem. Fundam.

    (1967)
  • D.F. James et al.

    The laminar flow of dilute polymer solutions through porous media

    J. Fluid Mech.

    (1975)
  • G.H. McKinley et al.

    Nonlinear dynamics of viscoelastic flow in axisymmetric abrupt contractions

    J. Fluid Mech.

    (1991)
  • P. Pakdel et al.

    Elastic instability and curved streamlines

    Phys. Rev. Lett.

    (1996)
  • Cited by (45)

    • Interface-resolved simulations of the confinement effect on the sedimentation of a sphere in yield-stress fluids

      2022, Journal of Non-Newtonian Fluid Mechanics
      Citation Excerpt :

      The dimensionless parameters employed for the simulations presented here are reported in Table 1. The present three-dimensional numerical solver has been utilized and extensively validated in the past for particulate flows [43–45], non-Newtonian flows [16,42,46–48] and multiphase problems in non-Newtonian fluids [49]. The code has also been recently validated for suspensions of rigid and soft particles and droplets in EVP and viscoelastic fluids [39].

    • Deformation of an axisymmetric viscoplastic drop in extensional/compressional flow

      2021, Journal of Non-Newtonian Fluid Mechanics
      Citation Excerpt :

      At relatively high deformation, shown in Fig. 8(b), the radial un-yielded torus shows diminished velocity and low shear intensity. Similar characteristics are discussed by De Vita et al. [42] and Chaparian et al. [43] studying visco-plastic flow in porous media. We conclude this section by providing phase planes that show the boundary between the region of stable stationary shapes and the region beyond the loss of stability.

    View all citing articles on Scopus
    View full text