General hydrodynamic features of elastoviscoplastic fluid flows through randomised porous media

Abstract A numerical study of yield-stress fluids flowing in porous media is presented. The porous media is randomly constructed by non-overlapping mono-dispersed circular obstacles. Two class of rheological models are investigated: elastoviscoplastic fluids (i.e. Saramito model) and viscoplastic fluids (i.e. Bingham model). A wide range of practical Weissenberg and Bingham numbers is studied at three different levels of porosities of the media. The emphasis is on revealing some physical transport mechanisms of yield-stress fluids in porous media when the elastic behaviour of this kind of fluids is incorporated. Thus, computations of elastoviscoplastic fluids are performed and are compared with the viscoplastic fluid flow properties. At a constant Weissenberg number, the pressure drop increases both with the Bingham number and the solid volume fraction of obstacles. However, the effect of elasticity is less trivial. At low Bingham numbers, the pressure drop of an elastoviscoplastic fluid increases compared to a viscoplastic fluid, while at high Bingham numbers we observe drag reduction by elasticity. At the yield limit (i.e. infinitely large Bingham numbers), elasticity of the fluid systematically promotes yielding: elastic stresses help the fluid to overcome the yield stress resistance at smaller pressure gradients. We observe that elastic effects increase with both Weissenberg and Bingham numbers. In both cases, elastic effects finally make the elastoviscoplastic flow unsteady, which consequently can result in chaos and turbulence. Graphical abstract


Introduction
Fluid flow through porous media is one of the classic fluid mechanics problems due to its importance for many industrial and environmental applications.Yield-stress fluids flowing in porous media is also intrinsic to a wide range of processes from cement grouting, enhanced oil recovery, hydraulic fracturing, etc. in construction and oil & gas industries to vertebroplasty (i.e.injection of bone "cement" into fractured vertebra) [1].
Over the past decades, major forward steps have been made in revealing the features of the flow by connecting the micro-scale features to the macro-scale hydrodynamic properties by integrating the linear equations governing the viscous fluids based on Darcy's findings.His framework was a baseline for further progressions in the fluid mechanics of porous media when the working fluid is Newtonian.
However, following experiments and theoretical works highlighted the invalidity of such approaches when it comes to more complex fluids such as viscoelastic and yieldstress fluids due to the non-linearity of the constitutive equations.More specifically for the yield-stress fluids, which are the focus of the present study, the issue is two-fold: (i) Yield-stress fluids rheologically behave like solids when the applied stress is less than a threshold-the yield stress.If we translate this to porous media terminology, it means that if the applied pressure gradient is less than a threshold, then there is no flow inside the porous media.Hence, a finite pressure gradient is required to initiate the flow of a yield-stress fluid inside a porous medium.(ii) Beyond the above limit (i.e. the yield limit), due to non-linearity of the constitutive equations, Darcy's approach is not valid.Hence, permeability does not only depend on the geometrical features of the porous media (e.g.porosity and shape of the obstacles/void spaces), but also it depends on the rheological parameters of the fluid.Consequently, we cannot decouple fluid's rheology and pore geometry which is feasible in Darcy's approach for Newtonian fluids.
These features have been demonstrated in several theoretical, numerical and experimental studies [2][3][4][5][6][7] and three main regimes have been identified for yield-stress fluid flows through porous media knowing that when the applied pressure gradient is less than the critical pressure gradient, there is no flow.If the applied pressure gradient slightly exceeds the critical value, the flow is extremely localised and goes through a single channel [5,8].In this regime, the flow rate linearly scales with the excessive pressure gradient.As the applied pressure gradient increases, a second regime emerges in which more and more channels will appear and the flow rate scales quadratically with the excessive applied pressure gradient.In the third regime, where the applied pressure gradient is much higher than the critical value, the flow rate again scales linearly with the excessive pressure gradient.Talon and co-workers in a series of studies [2,3,9] distinguished these three distinct regimes in the yield-stress fluid flow through porous media, which have been validated further by Chaparian & Tammisola [6] in randomised porous media.Emerging sophisticated experimental techniques in the recent years have revealed that "simple" yield-stress fluid rheological models are not sufficient for describing complex hydrodynamic features of many practical yield-stress fluids [10][11][12].Therefore, a number of rheological models have been proposed which include thixotropic and/or elastic behaviour of the fluid as well.A minimal elastoviscoplastic model is proposed by Saramito [13] which is indeed a combination of the Bingham and the Oldroyd-B models and is capable of capturing more crucial facets of complex yield-stress fluid flows [12,14,15].Hence, we use this rheological model in the present study.
We previously analysed the elastoviscoplastic fluid flow at micro-scale inside model porous media, i.e. flow over obstacles in a single periodic cell [4,16].We found that for relatively small Bingham numbers, the slightly elastic behaviour of the fluid increases the pressure gradient.However, at moderate Bingham numbers, the pressure gradients for viscoplastic and weakly elastoviscoplastic fluids were found to be very similar.The elasticity, nevertheless, changes the shape of unyielded and fouling regions to some extent.Indeed, in symmetric model geometries, given that the inertia is negligible, the unyielded/fouling regions shape is symmetric for viscoplastic fluids, while it is asymmetric for elastoviscoplastic fluids.This fact has been thoroughly discussed in [12,17].
Here, the focus is again on elastoviscoplastic fluid flows in porous media.The objectives are: (i) pushing the analysis further to a wider range of the Bingham numbers (i.e.yield limit) and porosities; (ii) addressing the effect of the Weissenberg number (i.e. the extent of elastic behaviour of the yield-stress fluid) and more importantly, (iii) investigating the randomised porous media (rather than a periodic pore-scale cell) and also bridging our previous micro-scale study to macro-scale features.Moreover, we shed some light on how the elastic effects make the flow oscillatory by time which can eventually result in elastic instabilities.
Elastic turbulence found some growing attention during the past decades [18,19] which is attributed to the influence of polymer chains in the creeping flow regime where the inertial forces are negligible.This results in chaotic patterns which are absent in corresponding creeping flow of Newtonian fluids.In elastic turbulence, flow behavior is mainly governed by non-linear elastic stresses characterised by large Weissenberg numbers (1 ≪ W i) or consequently large elasticity numbers which has a wide range of applications, spanning from materials science to the biomedical engineering [18,20].
For several yield-stress fluid flows, it has been reported that high elastic effects can be observed even at low Weissenberg numbers if the Bingham number is large enough [17,21,22].This fact resulted in proposing W i × Bn as the effective parameter controlling the emergence of elastic effects by Chaparian & Tammisola [17] which has been experimentally demonstrated as well by Villalba et al. [12].Indeed, synergy of elasticity and plasticity of the fluid triggers elastic effects even at small Weissenberg numbers.Hence, we expect that alternative "turbulence" happens in elastoviscoplastic fluid flows where both Reynolds and Weissenberg numbers are small while the Bingham number is fairly large.This has been slightly discussed here for flows in porous media although it is not the main aim of the present study.
In what follows, we set the problem in section 2 and briefly discuss the utilized numerical methods.The results are presented in section 3 and conclusions & discussions in section 4.

Problem setup
In the present study, we investigate two-dimensional very low Reynolds number flows of yield-stress fluids (viscoplastic & elastoviscoplastic) through porous media composed of mono-dispersed rigid circular obstacles (X) in a square domain (Ω) of size L × L. In this section, the governing equations, porous media construction and numerical methods are explained.To model the rheology of the elastoviscoplastic material, here we rely on Saramito model [23], in which the material is described as a Kelvin-Voigt viscoelastic solid before yielding, and a Bingham fluid with an extra elastic memory after yielding.

Mathematical formulation
The non-dimensional continuity, momentum and constitutive equations are as follows: where t is time, u = (u, v) is the velocity vector, p denotes the pressure, τ the stress tensor (sometimes is termed as extra stress tensor), and γ = (∇u + ∇u T ) is the rate of deformation tensor where ∇u = ∂u i /∂x j .Here, ∥ ⋅ ∥ v is the equivalent von Mises stress, (⋅) + is the positive part of the argument (i.e. is equal to zero if the argument is negative and otherwise is equal to the argument), and ▽ ⋅ is the upper-convected derivative: In these equations, the velocity vector is scaled with Û (i.e. the mean inlet velocity), pressure and stress tensor is scaled with the characteristic viscous stress μ Û / D where D (i.e.diameter of the obstacles) is the length scale and μ is the total viscosity.Moreover, β = μs /μ where μs is the solvent viscosity.
The non-dimensional parameters are the Reynolds number, Weissenberg number, and the Bingham number: The Reynolds number throughout this study is fixed and is equal to unity.The Weissenberg number represents as the ratio between elastic and viscous forces, and the Bingham number describes the ratio of the yield stress of the material to the characteristic viscous stress.In Appendix A, we discuss other possible scalings.Various approaches have been used in the numerical simulation of viscoelastic polymer solutions to handle the numerical instabilities appearing at "high" Weissenberg numbers [20,[24][25][26][27].In this study, the log-conformation approach, proposed by Fattal and Kupferman [28,29], is used.Although W i range studied here is not "high", in elastoviscoplastic fluids the elastic features can appear even at relatively low Weissenberg numbers, see [12,17,30].For the case of the Saramito model (i.e.equation 3), the conformation tensor A reads, For more details on log-conformation formulation, the readers are refereed to [31].

Porous medium construction
As mentioned, here we are interested in 2D flows through porous media constructed of mono-dispersed rigid circular obstacles (X) in a domain (Ω) of size L × L. These obstacles have a unit diameter (i.e. the length scale used in the above equations is the diameter of the obstacles) and L = 24.The centres of obstacles are chosen entirely random with uniform distribution.The constraints are non-overlapping and also being completely in the computational box (i.e. they do not touch the borders of the computational domain).
A sketch illustrating the considered geometries is presented in Fig. 1.The "volume fraction" of the obstacles is denoted by ϕ = meas(X)/meas(Ω) = N π/4L 2 where N is the number of obstacles.Hence, the porosity of the medium Φ (i.e.void fraction) can be defined simply as Φ = 1−ϕ.In the present study, we study three cases with the solid volume fraction 10%, 20%, and 30% (i.e.porosity Φ = 90, 80 and 70%), as displayed in Fig. 1.
Although in the present study, due to high computational cost (especially the elastoviscoplastic problem), we do not follow a statistical approach to simulate the flow in many realisations, the position of the obstacles are chosen completely random in the studied cases to avoid any bias in the results.Indeed, all the results are reported based on the realisations shown in Fig. 1 in which the obstacles are randomly generated to mimic the desired porosity.

Highlights of the computational details
In the following simulations, the no-slip boundary condition is imposed at the surface of the obstacles (u = 0 on ∂X), and free-slip boundary conditions at the lateral (i.e.top and bottom) boundaries, and the periodic boundary condition is implemented in the streamwise direction, while the left boundary is the inlet (i.e. the fluid flows from left to right).
In what follows, we briefly explain the two computational methods used in the present study to solve elastociscoplastic & viscoplastic problems.For more details, the readers are referred to our previous studies where we developed/implemented these methods.These references are mentioned accordingly.

Viscoplastic fluid flow simulations
We implement augmented Lagrangian method to simulate viscoplastic fluid flows [32,33].In this method, regularisation is avoided as the method is capable of handling the discontinuity in the Bingham model by relaxing the rate of the strain tensor.An open source finite element environment FreeFEM++ [34] is used for discretization and meshing which has been widely validated in our previous studies [6,[35][36][37][38].Anisotropic adaptive mesh in Ω/ X is combined with this method to get smoother yield surfaces and ensure high resolution of the flow features [17,32].

Elastoviscoplastic fluid flow simulations
The evolution of the conformation tensor (i.e.Eqs. 3 & 4) and the momentum equation are discritized using second-order central finite difference scheme on a uniform staggered grid in Ω, except for the advection term in Eq. 4 for which we use a fifth-order weighted essentially non-oscillatory (WENO) scheme [39,40].The time integration is performed with a fractional-step, third-order explicit Runge-Kutta scheme [41] for all equations.For more detail, the reader is referred to [22,42].
In the present study, the immersed boundary method (IBM) is used to represent the rigidity of obstacles (i.e.X), coupled with the uniform computational grid Ω, and hence exploit the possibility of utilizing efficient computational algorithms such as highly-scalable, fast Fourier transform (FFT)-based pressure solver.Hence in this way, the simulation is performed on a uniform Cartesian Eulerian grid with the solid surface reproduced by a uniformly distributed Lagrangian grid [42][43][44][45].The Eulerian grid is uniform in the two directions of the computational domain and the Lagrangian points are equally distributed on the solid surfaces-32 points per diameter equivalent to 101 points on the surface of each obstacle.The in-house developed IBM code has been widely verified in our previous physical studies involving multiphase elastoviscoplastic fluid flows such as flows in porous media [16,46], particle migration [30], and flow past an obstacle [42,45,47].
For the present flow configuration, we provide an estimate of the drag coefficient around a single cylinder in viscoelastic flow obtained with the present IBM method, for time-independent and time-dependent flow cases.More comprehensive results are available in [45].It should be mentioned that some previous studies indicate that the accuracy of the IBM method in viscoelastic flow can sometimes be improved by the so-called smooth continuation which allows a higher accuracy of stress gradients at the boundary [50,51], and at W i = 0.1 such method could provide an error of < 0.1%.However, the accuracy of the smooth IBM approach seemed to deteriorate with increasing Weissenberg number, and already at W i = 1.4 (W i = 0.7 by their definition) the error is 4.4% with 128 points per cylinder surface (see Fig. 7 in [51]), which is comparable to our IBM method error at W i = 2.In this study, we also refrain from the use of artificial diffusion, which may have significantly affected elastic instabilities [52] and would have been necessary if otherwise accurate spectral discretisations were employed as in [51].

Results
The emphasis in the present study is on studying the effect of the Weissenberg and Bingham numbers, and also the porosity of the medium on the hydrodynamic features.Hence, we fix β = 0.5 in all the elastoviscoplastic simulations and to minimize the inertial effects, we specify Re ≈ 1.Nevertheless, we vary the Bingham number in the range of 0 to 500.Knowing that the relaxation times of most practical yield-stress fluids are relatively small compared to the viscoelastic fluids [12,53], we simulate the flows at relatively small Weissenberg numbers, W i = 0.05 & 0.5, except at ϕ = 10% where the extreme case W i = 10 is also studied.The velocity contours (i.e.|u|) for the ϕ = 10%, 20%, 30% cases at W i = 0.5 and Bn = [0, 10, 500] are presented in Fig. 2. Each row represents a constant solid volume fraction (i.e.ϕ) and it increases from top to bottom.Also, the Bingham number increases from the left to the right columns while the Weissenberg number is fixed at 0.5.Fig. 3 shows again the velocity contour but the emphasis is on the effect of the Weissenberg number on the flow, hence at ϕ = 10%, a wide range of W i is covered for different Bingham numbers.
An increase in the Bingham number from Bn = 0 to Bn = 500 results in a more confined flow or "channelisation" which is well-documented before in the yield-stress fluids literature [2,5,6,8].Indeed, at high Bingham numbers, a large portion of the fluid remains unyielded and the flow emerges in thin channels through the porous media due to the high yield stress of the fluid compared to the deriving stress imposed by the applied pressure gradient.
Although changing the Weissenberg number at a fixed Bingham number does not have significant effect on the contour of velocity (see Fig. 3), this increase results in appearing high velocity regions in the open channels when the fluid passes through through the porous media due to the high yield stress of the fluid compared to the deriving stress imposed by the applied pressure gradient.
Although changing the Weissenberg number at a fixed Bingham number does not have significant effect on the contour of velocity (see Fig. 3), this increase results in appearing high velocity regions in the open 175 channels when the fluid passes through the gaps between the obstacles which is intuitive for more elastic fluids.This becomes clear by comparing the top and bottom rows where the yellow and the red contours appear in the middle of the open channels for high W i flows (please note that the same colour ranges are used for the top and bottom counterparts).
To further analyse the effects of the Weissenberg and Bingham numbers on the flow features, Figs.The trace of the conformation tensor(Fig.4) is the sign of the elongation of the polymer chains in viscoelastic fluids.Knowing that the trace is zero for Newtonian and viscoplastic fluids (i.e. the stress tensor is deviatoric), for elastoviscoplastic fluids, the trace of the stress tensor has the same interpretation as the viscoelastic fluids: the extent of stress generated by the elongation of fluid elements due to its elastic stress difference (i.e.N 1 = τ xx − τ yy ), respectively.The trace of the stress tensor (Fig. 4) is the sign of the elongation of the polymer chains in viscoelastic fluids.Knowing that the trace is zero for Newtonian and viscoplastic fluids (i.e. the stress tensor is deviatoric), for elastoviscoplastic fluids, the trace of the stress tensor has the same interpretation as the viscoelastic 185 fluids: the extent of stress generated by the elongation of fluid elements due to its elastic behaviour.As seen in Fig. 4, increasing the Weissenberg number from 0.05 (top panels) to 10 (bottom panels) results in an increase in the trace and also appearance of long and thin streaks in the streamwise direction.Indeed, longer streaks of high trace regions extend further downstream from the cylinders' surface.Moreover, increasing the Bingham number from Bn = 0 to Bn = 500 while the Weissenberg number is kept constant (compare 190 top or bottom panels from left to right) leads to a more pronounced extension of fluid particles as indicated by the increase in the maximum of the trace, both locally and in the entire domain.For instance, when W i = 0.05 (top panels of Fig. 4), increasing the Bingham number from 0 to 500 changes the trace by almost three orders of magnitude: from dark blue (left panel) to yellow (right panel) contours.Notably, keeping the Weissenberg number the same and increasing the Bingham number leads to more elastic responses which is 195 previously reported in [11,16] for complex internal flows of elastoviscoplastic fluids.It should be noted that also increasing the porosity leads to stronger streaks which is not shown here for the sake of brevity.
As mentioned, Fig. 5 depicts the contours of N 1 .In the current study, it can be interpreted as a sign of localisation since indeed τ xx is much higher than τ yy when N 1 is large.As it is clear from Fig. 5, when B = 0 (the Newtonian case), the entire fluid flows and N 1 is approximately negligible.However, at the other 200 extreme, when B = 500 (i.e.moving towards the yield limit) and the flow is highly localised in one or few open channels, N 1 is quite large.It is also shown that the magnitude of N 1 increases with the increasing solid volume fraction.The value of N 1 is positive throughout the domain, except in the vicinity of the cylinders' upstream stagnation points.This can be interpreted as the transition of the flow direction from 7 behaviour.As seen in Fig. 4, increasing the Weissenberg number from 0.05 (top panels) to 10 (bottom panels) results in an increase in the trace and also appearance of long and thin streaks in the streamwise direction.Indeed, longer streaks of high trace regions extend further downstream from the cylinders' surface.Moreover, increasing the Bingham number from Bn = 0 to Bn = 500 while the Weissenberg number is kept constant (compare top or bottom panels from left to right) leads to a more pronounced extension of fluid particles as indicated by the increase in the maximum of the trace, both locally and in the entire domain.For instance, when W i = 0.05 (top panels of Fig. 4), increasing the Bingham number from 0 to 500 changes the trace by almost three orders of magnitude: from dark blue (left panel) to yellow (right panel) contours.Notably, keeping the Weissenberg number the same and increasing the Bingham number leads to more elastic responses which is previously reported in [12,17] for complex internal flows of elastoviscoplastic fluids.It should be noted that also increasing the porosity leads to stronger streaks which is not shown here for the sake of brevity.
As mentioned, Fig. 5 depicts the contours of τ xx − τ yy .In the current study, it can be interpreted as a sign of localisation since indeed τ xx is much higher than τ yy when the normal stress difference is large.As it is clear from Fig. 5, when Bn = 0 (the Newtonian case), the entire fluid flows and τ xx − τ yy is approximately negligible.However, at the other extreme, when Bn = 500 (i.e.moving towards the yield limit) and the flow is highly localised in one or few open channels, τ xx − τ yy is quite large.It is also shown that the magnitude of τ xx − τ yy increases with the increasing solid To help us analyse the bulk transport properties, Fig. 6 shows the pressure gradient as a function of the Bingham number for di↵erent porosities (with di↵erent colours) and Weissenberg numbers (with di↵erent symbols & line styles; see the caption).First of all, it is clear that for higher volume fractions the pressure 210 drop is higher at the same W i and Bn numbers, which is intuitive since the fluid should be percolated in a more packed arrangement of obstacles.At small Bingham numbers, the least studied Weissenberg number flow (W i = 0.05) and the viscoplastic flow are indistinguishable.However, as we increase the Bingham number and keep the Weissenberg number fixed at 0.05, we see that the discrepancy increases which has been evidenced in a number of our previous studies: plasticity of the fluid triggers the elastic behaviour of 215 the material or indeed as W i × Bn increases, we observe more elastic features in the flow [11,16].
At lower Bn values, if we increase the Weissenberg number, the pressure gradient increases; this is in agreement with our previous study in a small periodic cell, i.e. [3].At low Bingham numbers, the flow is rather evenly distributed over the porous medium, i.e. the flow is not heterogeneous/localised and there are less unyielded/fouling regions.Hence, the drag increase as a function of the elasticity of the fluid at 220 low Bingham numbers can be attributed to the non-straight streamlines where more energy is required to percolate an elastic fluid compared to a non-elastic one (due to stored elastic energy in fluid elements; see [53]).On the other hand, if we now focus on the high Bingham numbers regime, the e↵ect of the elasticity of the fluid is inverse: indeed at high Bingham number regime, higher Weissenberg flows have less pressure gradient compared to the small Weissenberg flows.This is the regime which was not explored in our previous 225 study [3].The reason for drag reduction with elasticity at high Bingham numbers can be attributed to the fact that at higher Bingham numbers, the flow becomes more and more channelised (i.e.localised in one 8 To help us analyse the bulk transport properties, Fig. 6 shows the pressure gradient as a function of the Bingham number for di↵erent porosities (with di↵erent colours) and Weissenberg numbers (with di↵erent symbols & line styles; see the caption).First of all, it is clear that for higher volume fractions the pressure 210 drop is higher at the same W i and Bn numbers, which is intuitive since the fluid should be percolated in a more packed arrangement of obstacles.At small Bingham numbers, the least studied Weissenberg number flow (W i = 0.05) and the viscoplastic flow are indistinguishable.However, as we increase the Bingham number and keep the Weissenberg number fixed at 0.05, we see that the discrepancy increases which has been evidenced in a number of our previous studies: plasticity of the fluid triggers the elastic behaviour of 215 the material or indeed as W i × Bn increases, we observe more elastic features in the flow [11,16].
At lower Bn values, if we increase the Weissenberg number, the pressure gradient increases; this is in agreement with our previous study in a small periodic cell, i.e. [3].At low Bingham numbers, the flow is rather evenly distributed over the porous medium, i.e. the flow is not heterogeneous/localised and there are less unyielded/fouling regions.Hence, the drag increase as a function of the elasticity of the fluid at 220 low Bingham numbers can be attributed to the non-straight streamlines where more energy is required to percolate an elastic fluid compared to a non-elastic one (due to stored elastic energy in fluid elements; see [53]).On the other hand, if we now focus on the high Bingham numbers regime, the e↵ect of the elasticity of the fluid is inverse: indeed at high Bingham number regime, higher Weissenberg flows have less pressure gradient compared to the small Weissenberg flows.This is the regime which was not explored in our previous 225 study [3].The reason for drag reduction with elasticity at high Bingham numbers can be attributed to the fact that at higher Bingham numbers, the flow becomes more and more channelised (i.e.localised in one 8 To help us analyse the bulk transport properties, Fig. 6 shows the pressure gradient as a function of the Bingham number for di↵erent porosities (with di↵erent colours) and Weissenberg numbers (with di↵erent symbols & line styles; see the caption).First of all, it is clear that for higher volume fractions the pressure 210 drop is higher at the same W i and Bn numbers, which is intuitive since the fluid should be percolated in a more packed arrangement of obstacles.At small Bingham numbers, the least studied Weissenberg number flow (W i = 0.05) and the viscoplastic flow are indistinguishable.However, as we increase the Bingham number and keep the Weissenberg number fixed at 0.05, we see that the discrepancy increases which has been evidenced in a number of our previous studies: plasticity of the fluid triggers the elastic behaviour of 215 the material or indeed as W i × Bn increases, we observe more elastic features in the flow [11,16].
At lower Bn values, if we increase the Weissenberg number, the pressure gradient increases; this is in agreement with our previous study in a small periodic cell, i.e. [3].At low Bingham numbers, the flow is rather evenly distributed over the porous medium, i.e. the flow is not heterogeneous/localised and there are less unyielded/fouling regions.Hence, the drag increase as a function of the elasticity of the fluid at 220 low Bingham numbers can be attributed to the non-straight streamlines where more energy is required to percolate an elastic fluid compared to a non-elastic one (due to stored elastic energy in fluid elements; see [53]).On the other hand, if we now focus on the high Bingham numbers regime, the e↵ect of the elasticity of the fluid is inverse: indeed at high Bingham number regime, higher Weissenberg flows have less pressure gradient compared to the small Weissenberg flows.This is the regime which was not explored in our previous 225 study [3].The reason for drag reduction with elasticity at high Bingham numbers can be attributed to the fact that at higher Bingham numbers, the flow becomes more and more channelised (i.e.localised in one 8 To help us analyse the bulk transport properties, Fig. 6 shows the pressure gradient as a function of the Bingham number for di↵erent porosities (with di↵erent colours) and Weissenberg numbers (with di↵erent symbols & line styles; see the caption).First of all, it is clear that for higher volume fractions the pressure 210 drop is higher at the same W i and Bn numbers, which is intuitive since the fluid should be percolated in a more packed arrangement of obstacles.At small Bingham numbers, the least studied Weissenberg number flow (W i = 0.05) and the viscoplastic flow are indistinguishable.However, as we increase the Bingham number and keep the Weissenberg number fixed at 0.05, we see that the discrepancy increases which has been evidenced in a number of our previous studies: plasticity of the fluid triggers the elastic behaviour of 215 the material or indeed as W i × Bn increases, we observe more elastic features in the flow [11,16].
At lower Bn values, if we increase the Weissenberg number, the pressure gradient increases; this is in agreement with our previous study in a small periodic cell, i.e. [3].At low Bingham numbers, the flow is rather evenly distributed over the porous medium, i.e. the flow is not heterogeneous/localised and there are less unyielded/fouling regions.Hence, the drag increase as a function of the elasticity of the fluid at 220 low Bingham numbers can be attributed to the non-straight streamlines where more energy is required to percolate an elastic fluid compared to a non-elastic one (due to stored elastic energy in fluid elements; see [53]).On the other hand, if we now focus on the high Bingham numbers regime, the e↵ect of the elasticity of the fluid is inverse: indeed at high Bingham number regime, higher Weissenberg flows have less pressure gradient compared to the small Weissenberg flows.This is the regime which was not explored in our previous 225 study [3].The reason for drag reduction with elasticity at high Bingham numbers can be attributed to the fact that at higher Bingham numbers, the flow becomes more and more channelised (i.e.localised in one 8 To help us analyse the bulk transport properties, Fig. 6 shows the pressure gradient as a function of the Bingham number for di↵erent porosities (with di↵erent colours) and Weissenberg numbers (with di↵erent symbols & line styles; see the caption).First of all, it is clear that for higher volume fractions the pressure 210 drop is higher at the same W i and Bn numbers, which is intuitive since the fluid should be percolated in a more packed arrangement of obstacles.At small Bingham numbers, the least studied Weissenberg number flow (W i = 0.05) and the viscoplastic flow are indistinguishable.However, as we increase the Bingham number and keep the Weissenberg number fixed at 0.05, we see that the discrepancy increases which has been evidenced in a number of our previous studies: plasticity of the fluid triggers the elastic behaviour of 215 the material or indeed as W i × Bn increases, we observe more elastic features in the flow [11,16].
At lower Bn values, if we increase the Weissenberg number, the pressure gradient increases; this is in agreement with our previous study in a small periodic cell, i.e. [3].At low Bingham numbers, the flow is rather evenly distributed over the porous medium, i.e. the flow is not heterogeneous/localised and there are less unyielded/fouling regions.Hence, the drag increase as a function of the elasticity of the fluid at 220 low Bingham numbers can be attributed to the non-straight streamlines where more energy is required to percolate an elastic fluid compared to a non-elastic one (due to stored elastic energy in fluid elements; see [53]).On the other hand, if we now focus on the high Bingham numbers regime, the e↵ect of the elasticity of the fluid is inverse: indeed at high Bingham number regime, higher Weissenberg flows have less pressure gradient compared to the small Weissenberg flows.This is the regime which was not explored in our previous 225 study [3].The reason for drag reduction with elasticity at high Bingham numbers can be attributed to the fact that at higher Bingham numbers, the flow becomes more and more channelised (i.e.localised in one 8 streamwise to vertical direction in order to align with the surface of the solid cylinders.Larger Bingham 205 numbers and solid volume fractions increase the stagnation points e↵ects and leads to an increase in the area with N 1 0 shown by the blue colour contours. To help us analyse the bulk transport properties, Fig. 6 shows the pressure gradient as a function of the Bingham number for di↵erent porosities (with di↵erent colours) and Weissenberg numbers (with di↵erent symbols & line styles; see the caption).First of all, it is clear that for higher volume fractions the pressure 210 drop is higher at the same W i and Bn numbers, which is intuitive since the fluid should be percolated in a more packed arrangement of obstacles.At small Bingham numbers, the least studied Weissenberg number flow (W i = 0.05) and the viscoplastic flow are indistinguishable.However, as we increase the Bingham number and keep the Weissenberg number fixed at 0.05, we see that the discrepancy increases which has been evidenced in a number of our previous studies: plasticity of the fluid triggers the elastic behaviour of 215 the material or indeed as W i × Bn increases, we observe more elastic features in the flow [11,16].
At lower Bn values, if we increase the Weissenberg number, the pressure gradient increases; this is in agreement with our previous study in a small periodic cell, i.e. [3].At low Bingham numbers, the flow is rather evenly distributed over the porous medium, i.e. the flow is not heterogeneous/localised and there are less unyielded/fouling regions.Hence, the drag increase as a function of the elasticity of the fluid at 220 low Bingham numbers can be attributed to the non-straight streamlines where more energy is required to percolate an elastic fluid compared to a non-elastic one (due to stored elastic energy in fluid elements; see [53]).On the other hand, if we now focus on the high Bingham numbers regime, the e↵ect of the elasticity of the fluid is inverse: indeed at high Bingham number regime, higher Weissenberg flows have less pressure gradient compared to the small Weissenberg flows.This is the regime which was not explored in our previous 225 study [3].The reason for drag reduction with elasticity at high Bingham numbers can be attributed to the fact that at higher Bingham numbers, the flow becomes more and more channelised (i.e.localised in one 8 streamwise to vertical direction in order to align with the surface of the solid cylinders.Larger Bingham 205 numbers and solid volume fractions increase the stagnation points e↵ects and leads to an increase in the area with N 1 0 shown by the blue colour contours.To help us analyse the bulk transport properties, Fig. 6 shows the pressure gradient as a function of the Bingham number for di↵erent porosities (with di↵erent colours) and Weissenberg numbers (with di↵erent symbols & line styles; see the caption).First of all, it is clear that for higher volume fractions the pressure 210 drop is higher at the same W i and Bn numbers, which is intuitive since the fluid should be percolated in a more packed arrangement of obstacles.At small Bingham numbers, the least studied Weissenberg number flow (W i = 0.05) and the viscoplastic flow are indistinguishable.However, as we increase the Bingham number and keep the Weissenberg number fixed at 0.05, we see that the discrepancy increases which has been evidenced in a number of our previous studies: plasticity of the fluid triggers the elastic behaviour of 215 the material or indeed as W i × Bn increases, we observe more elastic features in the flow [11,16].
At lower Bn values, if we increase the Weissenberg number, the pressure gradient increases; this is in agreement with our previous study in a small periodic cell, i.e. [3].At low Bingham numbers, the flow is rather evenly distributed over the porous medium, i.e. the flow is not heterogeneous/localised and there are less unyielded/fouling regions.Hence, the drag increase as a function of the elasticity of the fluid at 220 low Bingham numbers can be attributed to the non-straight streamlines where more energy is required to percolate an elastic fluid compared to a non-elastic one (due to stored elastic energy in fluid elements; see [53]).On the other hand, if we now focus on the high Bingham numbers regime, the e↵ect of the elasticity of the fluid is inverse: indeed at high Bingham number regime, higher Weissenberg flows have less pressure gradient compared to the small Weissenberg flows.This is the regime which was not explored in our previous 225 study [3].The reason for drag reduction with elasticity at high Bingham numbers can be attributed to the fact that at higher Bingham numbers, the flow becomes more and more channelised (i.e.localised in one 8 streamwise to vertical direction in order to align with the surface of the solid cylinders.Larger Bingham 205 numbers and solid volume fractions increase the stagnation points e↵ects and leads to an increase in the area with N 1 0 shown by the blue colour contours. To help us analyse the bulk transport properties, Fig. 6 shows the pressure gradient as a function of the Bingham number for di↵erent porosities (with di↵erent colours) and Weissenberg numbers (with di↵erent symbols & line styles; see the caption).First of all, it is clear that for higher volume fractions the pressure 210 drop is higher at the same W i and Bn numbers, which is intuitive since the fluid should be percolated in a more packed arrangement of obstacles.At small Bingham numbers, the least studied Weissenberg number flow (W i = 0.05) and the viscoplastic flow are indistinguishable.However, as we increase the Bingham number and keep the Weissenberg number fixed at 0.05, we see that the discrepancy increases which has been evidenced in a number of our previous studies: plasticity of the fluid triggers the elastic behaviour of 215 the material or indeed as W i × Bn increases, we observe more elastic features in the flow [11,16].
At lower Bn values, if we increase the Weissenberg number, the pressure gradient increases; this is in agreement with our previous study in a small periodic cell, i.e. [3].At low Bingham numbers, the flow is rather evenly distributed over the porous medium, i.e. the flow is not heterogeneous/localised and there are less unyielded/fouling regions.Hence, the drag increase as a function of the elasticity of the fluid at 220 low Bingham numbers can be attributed to the non-straight streamlines where more energy is required to percolate an elastic fluid compared to a non-elastic one (due to stored elastic energy in fluid elements; see [53]).On the other hand, if we now focus on the high Bingham numbers regime, the e↵ect of the elasticity of the fluid is inverse: indeed at high Bingham number regime, higher Weissenberg flows have less pressure gradient compared to the small Weissenberg flows.This is the regime which was not explored in our previous 225 study [3].The reason for drag reduction with elasticity at high Bingham numbers can be attributed to the fact that at higher Bingham numbers, the flow becomes more and more channelised (i.e.localised in one 8 streamwise to vertical direction in order to align with the surface of the solid cylinders.Larger Bingham 205 numbers and solid volume fractions increase the stagnation points e↵ects and leads to an increase in the area with N 1 0 shown by the blue colour contours. To help us analyse the bulk transport properties, Fig. 6 shows the pressure gradient as a function of the Bingham number for di↵erent porosities (with di↵erent colours) and Weissenberg numbers (with di↵erent symbols & line styles; see the caption).First of all, it is clear that for higher volume fractions the pressure 210 drop is higher at the same W i and Bn numbers, which is intuitive since the fluid should be percolated in a more packed arrangement of obstacles.At small Bingham numbers, the least studied Weissenberg number flow (W i = 0.05) and the viscoplastic flow are indistinguishable.However, as we increase the Bingham number and keep the Weissenberg number fixed at 0.05, we see that the discrepancy increases which has been evidenced in a number of our previous studies: plasticity of the fluid triggers the elastic behaviour of 215 the material or indeed as W i × Bn increases, we observe more elastic features in the flow [11,16].
At lower Bn values, if we increase the Weissenberg number, the pressure gradient increases; this is in agreement with our previous study in a small periodic cell, i.e. [3].At low Bingham numbers, the flow is rather evenly distributed over the porous medium, i.e. the flow is not heterogeneous/localised and there are less unyielded/fouling regions.Hence, the drag increase as a function of the elasticity of the fluid at 220 low Bingham numbers can be attributed to the non-straight streamlines where more energy is required to percolate an elastic fluid compared to a non-elastic one (due to stored elastic energy in fluid elements; see [53]).On the other hand, if we now focus on the high Bingham numbers regime, the e↵ect of the elasticity of the fluid is inverse: indeed at high Bingham number regime, higher Weissenberg flows have less pressure gradient compared to the small Weissenberg flows.This is the regime which was not explored in our previous 225 study [3].The reason for drag reduction with elasticity at high Bingham numbers can be attributed to the fact that at higher Bingham numbers, the flow becomes more and more channelised (i.e.localised in one 8 streamwise to vertical direction in order to align with the surface of the solid cylinders.Larger Bingham 205 numbers and solid volume fractions increase the stagnation points e↵ects and leads to an increase in the area with N 1 0 shown by the blue colour contours.To help us analyse the bulk transport properties, Fig. 6 shows the pressure gradient as a function of the Bingham number for di↵erent porosities (with di↵erent colours) and Weissenberg numbers (with di↵erent symbols & line styles; see the caption).First of all, it is clear that for higher volume fractions the pressure 210 drop is higher at the same W i and Bn numbers, which is intuitive since the fluid should be percolated in a more packed arrangement of obstacles.At small Bingham numbers, the least studied Weissenberg number flow (W i = 0.05) and the viscoplastic flow are indistinguishable.However, as we increase the Bingham number and keep the Weissenberg number fixed at 0.05, we see that the discrepancy increases which has been evidenced in a number of our previous studies: plasticity of the fluid triggers the elastic behaviour of 215 the material or indeed as W i × Bn increases, we observe more elastic features in the flow [11,16].
At lower Bn values, if we increase the Weissenberg number, the pressure gradient increases; this is in agreement with our previous study in a small periodic cell, i.e. [3].At low Bingham numbers, the flow is rather evenly distributed over the porous medium, i.e. the flow is not heterogeneous/localised and there are less unyielded/fouling regions.Hence, the drag increase as a function of the elasticity of the fluid at 220 low Bingham numbers can be attributed to the non-straight streamlines where more energy is required to percolate an elastic fluid compared to a non-elastic one (due to stored elastic energy in fluid elements; see [53]).On the other hand, if we now focus on the high Bingham numbers regime, the e↵ect of the elasticity of the fluid is inverse: indeed at high Bingham number regime, higher Weissenberg flows have less pressure gradient compared to the small Weissenberg flows.This is the regime which was not explored in our previous 225 study [3].The reason for drag reduction with elasticity at high Bingham numbers can be attributed to the fact that at higher Bingham numbers, the flow becomes more and more channelised (i.e.localised in one 8 streamwise to vertical direction in order to align with the surface of the solid cylinders.Larger Bingham 205 numbers and solid volume fractions increase the stagnation points e↵ects and leads to an increase in the area with N 1 0 shown by the blue colour contours.To help us analyse the bulk transport properties, Fig. 6 shows the pressure gradient as a function of the Bingham number for di↵erent porosities (with di↵erent colours) and Weissenberg numbers (with di↵erent symbols & line styles; see the caption).First of all, it is clear that for higher volume fractions the pressure 210 drop is higher at the same W i and Bn numbers, which is intuitive since the fluid should be percolated in a more packed arrangement of obstacles.At small Bingham numbers, the least studied Weissenberg number flow (W i = 0.05) and the viscoplastic flow are indistinguishable.However, as we increase the Bingham number and keep the Weissenberg number fixed at 0.05, we see that the discrepancy increases which has been evidenced in a number of our previous studies: plasticity of the fluid triggers the elastic behaviour of 215 the material or indeed as W i × Bn increases, we observe more elastic features in the flow [11,16].
At lower Bn values, if we increase the Weissenberg number, the pressure gradient increases; this is in agreement with our previous study in a small periodic cell, i.e. [3].At low Bingham numbers, the flow is rather evenly distributed over the porous medium, i.e. the flow is not heterogeneous/localised and there are less unyielded/fouling regions.Hence, the drag increase as a function of the elasticity of the fluid at 220 low Bingham numbers can be attributed to the non-straight streamlines where more energy is required to percolate an elastic fluid compared to a non-elastic one (due to stored elastic energy in fluid elements; see [53]).On the other hand, if we now focus on the high Bingham numbers regime, the e↵ect of the elasticity of the fluid is inverse: indeed at high Bingham number regime, higher Weissenberg flows have less pressure gradient compared to the small Weissenberg flows.This is the regime which was not explored in our previous 225 study [3].The reason for drag reduction with elasticity at high Bingham numbers can be attributed to the fact that at higher Bingham numbers, the flow becomes more and more channelised (i.e.localised in one 8 volume fraction.The value of τ xx −τ yy is positive throughout the domain, except in the vicinity of the cylinders' upstream stagnation points.This can be interpreted as the transition of the flow direction from streamwise to vertical direction in order to align with the surface of the solid cylinders.Larger Bingham numbers and solid volume fractions increase the points effects and leads to an increase in the area with τ xx − τ yy ⩽ 0 shown by the blue colour contours. To help us analyse the bulk transport properties, Fig. 6 shows the pressure gradient as a function of the Bingham number for different porosities (with different colours) and Weissenberg numbers (with different symbols & line styles; see the caption).First of all, it is clear that for higher volume fractions the pressure drop is higher at the same W i and Bn numbers, which is intuitive since the fluid should be percolated in a more packed arrangement of obstacles.At small Bingham numbers, the least studied Weissenberg number flow (W i = 0.05) and the viscoplastic flow are indistinguishable.However, as we increase the Bingham number and keep the Weissenberg number fixed at 0.05, we see that the discrepancy increases which has been evidenced in a number of our previous studies: plasticity of the fluid triggers the elastic behaviour of the material or indeed as W i × Bn increases, we observe more elastic features in the flow [12,17].
At lower Bn values, if we increase the Weissenberg number, the pressure gradient increases; this is in agreement with our previous study in a small periodic cell, i.e. [4].At low Bingham numbers, the flow is rather evenly distributed over the porous medium, i.e. the flow is not heterogeneous/localised and there are less unyielded/fouling regions.Hence, the drag increase as a function of the elasticity of the fluid at low or few open channels) and hence the streamlines are more straight.Thus, a qualitative argument can be made based on a flow through a channel, even though the analogy is not perfect.In a laminar channel flow, the Fanning friction factor decreases with the Weissenberg number [21] (see [29] for the exact solution of Bingham numbers can be attributed to the non-straight streamlines where more energy is required to percolate an elastic fluid compared to a non-elastic one (due to stored elastic energy in fluid elements; see [54]).On the other hand, if we now focus on the high Bingham numbers regime, the effect of the elasticity of the fluid is inverse: indeed at high Bingham number regime, higher Weissenberg flows have less pressure gradient compared to the small Weissenberg flows.This is the regime which was not explored in our previous study [4].The reason for drag reduction with elasticity at high Bingham numbers can be attributed to the fact that at higher Bingham numbers, the flow becomes more and more channelised (i.e.localised in one or few open channels) and hence the streamlines are more straight.Thus, a qualitative argument can be made based on a flow through a channel, even though the analogy is not perfect.In a laminar channel flow, the Fanning friction factor decreases with the Weissenberg number [22] (see [30] for the exact solution of elastoviscoplastic fluid Poiseuille flow).
The main source of high drag in channel flow of yield-stress fluids is the sharp shear rate close to the walls which is mitigated by the elasticity of the fluid as the core plug region narrows down compared to the viscoplastic fluids.Thus, the elasticity of the fluid reduces the viscous stresses near the walls which results in less drag.More or less, the same mechanism can be observed here.In summary, we observed two different regimes for drag.There is a zone of Bingham number in which a transition from drag enhancement to drag reduction happens, for example, for the case ϕ = 10%, all curves associated with different W i meet around Bn ≈ 5. Obviously, this zone shifts towards higher Bingham numbers as the solid volume fraction is increasing: for ϕ = 20%, it is Bn ≈ 50 and for ϕ = 30% it is Bn ≈ 500.Nevertheless, it is clear that Y c is larger for higher Weissenberg flows.In other words, for the viscoplastic fluid (i.e.case W i = 0), ∆P /L scales with the Bingham number (see [4] for the proof) while the slope is much more gentle at higher Weissenberg numbers as Bn → ∞ (i.e. at the yield limit).This has been generally observed before in other physical problems such as particle sedimentation in elastoviscoplastic fluids [14]: elastic stresses aid the particle to overcome the yield stress and sediment.It should be noted that the data represented in Fig. 6 associated with elastoviscoplastic fluid is the time-averaged pressure gradient.The instantaneous pressure gradient (∆P /L versus time t) for different cases are plotted in Figs.7 and 8. Clearly, when the Weissenberg and the Bingham numbers are relatively small (left and middle panels in Fig. 7), the steady state solutions are achieved after rather short time intervals irrespective of the solid volume fraction.The right panels (c,f) in Fig. 7, however, show the cases of high Bingham number Bn = 500 (i.e.large W i × Bn).In these cases, the pressure gradient oscillates around an average value even after long times.Hence, in these cases, only a quasi-steady state is reached by the elastoviscoplastic fluid.The average values in the final quasi-steady state are previously reported in Fig. 6.As can be seen by comparing the different curves in Fig. 7 (panel (f) in particular), the oscillations' amplitude is larger when the solid volume fraction is higher.Nevertheless, the oscillation amplitudes compared to the average values of ∆P /L are rather negligible.Another point is that, in these cases, the flow converges to the quasi-steady state solutions faster when the solid volume fraction is smaller: for the case W i = 0.5 & Bn = 500 (see Fig. 7(f)), the flow reaches the quasi-steady state at t ≈ 5 when ϕ = 10%, at t = 10 when ϕ = 20% and at t = 20 when ϕ = 30%.
If we increase the value of W i × Bn even more (for instance W i = 10 & Bn = 500 which is represented in Fig. 8), there are cases in which the steady state solutions are not achievable, at least at practical time scales.This is especially pronounced at higher solid volume fractions.For example, the case W i = 10, Bn = 500 & ϕ = 0.2 which is marked in Fig. 8 by the red diamonds has not reached a plateau even at t = 200 and the slope also does not seem to decrease over the next time decades.This is the main reason that 10 ⩽ W i simulations are not reported in the present study when 0.1 < ϕ.

Conclusions and discussion
Numerical simulations of yield-stress fluid flows in porous media were performed.The focus was on comparing elastoviscoplastic (i.e.Saramito model [12]) and viscoplastic (i.e.Bingham fluid) rheological models.The porous media were constructed by randomised mono-dispersed non-overlapping circular obstacles where the size of the square computational domain is 48 times of the obstacle radius.Three different solid "volume" fractions (ϕ) are evaluated; ϕ = 0.1, 0.2 and 0.3.We presented a detailed analysis of the pressure drop for all of the case studies as a function of the Weissenberg number (W i), Bingham number (Bn) and the porosity of the media (1 − ϕ).
It has been demonstrated that ∆P /L increases with the Bingham number regardless of the Weissenberg and ϕ which is completely intuitive since the yield stress resistance increases.Also, the pressure drop increases with the solid volume fraction of the medium at fixed any W i and Bn.However, pressure drop has a non-monotonic behaviour with the Weissenberg number: while at the small Bingham numbers, ∆P /L increases with the Weissenberg number, at large Bn, the Weissenberg number has an opposite effect.In other words, at low Bingham numbers, elasticity increases the pressure drop, while at high Bingham numbers, elasticity has a drag reduction effect.This is systematically observed for all the considered solid volume fractions.
At the yield limit (Bn → ∞), the flow is highly localised to some channels while other parts of the fluid are quiescent.At this limit, ∆P /L ∼ Bn n where n is equal to one for "simple" viscoplastic fluids [3].However, as depicted in Fig. 6, for elastoviscoplastic fluids, n is less than one which renders the critical pressure gradient (or the inverse of the critical yield number; see [5,6]) to be smaller in elastoviscoplastic fluids compared to the viscoplastic fluids.
As previously observed in numerous number of studies on viscoelastic liquids, high elasticity of the fluid (i.e.high W i) results in unsteadiness, chaos and finally elastic turbulence in the inertialess regime [17]; see [54] in the context of porous media flows.In elastoviscoplastic fluid flows in a periodic cell (i.e.model porous media), we previously observed a transition to time-dependent oscillations at Bi = 10 and W i ≈ 0.5 [15], using a volume-penalization IBM approach with > 700 grid points over each cylinder surface.In  elastoviscoplastic fluids, the elastic behaviour is characterised by W i×Bn as evidenced both computationally and experimentally in our previous studies [11,16,20,21].Hence, even at relatively small Weissenberg numbers, if the Bingham number is sufficiently high, the unsteadiness appears; here we quantified this behaviour in Figs.7 and 8 where the pressure drop is shown as a function of time.Also, the effect of porosity has been discussed: at low solid volume fractions (e.g.ϕ = 10%), the flow is steady at low Bn and W i, while increasing the volume fraction to ϕ = 20% or 30% makes the flow unsteady at the same low 295 W i and Bn.In summary, the pressure drop is smaller and time-independent (steady) when the medium is more porous and the flow is at low Bn and W i numbers.Conversely, for larger values of the Bingham number, Weissenberg number and the solid volume fraction, the pressure drop becomes considerably timedependent.This is in line with studies of viscoelastic fluids in porous media, see e.g.[55].Recent studies [56] have suggested that the number of stagnation points in the flow affects the onset of elastic turbulence.

300
The number of stagnation points naturally increases with the solid volume fraction, but also depends on the precise geometric configuration in the disordered porous medium [56,57].Elastic turbulence results in a sudden increase in drag above a certain Weissenberg number threshold in viscoelastic flows through porous media [58].Hypothetically, if this were to happen in elastoviscoplastic flows as well when W i × Bn becomes very large, then our conclusion of drag reduction by elasticity may change.Characterising the elastoviscoplastic turbulence in porous media is, however, beyond the scope of the present work.Due to the high cost of numerical simulations of elastoviscoplastic fluid flows, here we focused on unravelling general features of this kind of fluids in porous media flows.Future studies are required to delve into the details and also the bulk transport properties which hopefully can yield general Darcy law for such kind of fluids.Thus, we only scratched the surface here which can be a baseline for future studies to extend our work in various directions, for example, focusing on more realistic topologies (e.g.digital rocks), more sophisticated rheological models (e.g.thixo-elasto-visco-plastic behaviours), or inertial regimes.Also, in industrial applications, the porosity of the medium can be very small such as "digital rocks" but for 2D randomised mono-dispersed obstacles, we cannot reach that limit.This remains for future studies to be addressed.

Conclusions and discussion
Numerical simulations of yield-stress fluid flows in porous media were performed.The focus was on comparing elastoviscoplastic (i.e.Saramito model [13]) and viscoplastic (i.e.Bingham fluid) rheological models.The porous media were constructed by randomised mono-dispersed non-overlapping circular obstacles where the size of the square computational domain is 48 times of the obstacle radius.Three different solid "volume" fractions (ϕ) are evaluated; ϕ = 0.1, 0.2 and 0.3.We presented a detailed analysis of the pressure drop for all of the case studies as a function of the Weissenberg number (W i), Bingham number (Bn) and the porosity of the media (1 − ϕ).
It has been demonstrated that ∆P /L increases with the Bingham number regardless of the Weissenberg and ϕ which is completely intuitive since the yield stress resistance increases.Also, the pressure drop increases with the solid volume fraction of the medium at fixed any W i and Bn.However, pressure drop has a non-monotonic behaviour with the Weissenberg number: while at the small Bingham numbers, ∆P /L increases with the Weissenberg number, at large Bn, the Weissenberg number has an opposite effect.In other words, at low Bingham numbers, elasticity increases the pressure drop, while at high Bingham numbers, elasticity has a drag reduction effect.This is systematically observed for all the considered solid volume fractions.
At the yield limit (Bn → ∞), the flow is highly localised to some channels while other parts of the fluid are quiescent.At this limit, ∆P /L ∼ Bn n where n is equal to one for "simple" viscoplastic fluids [4].However, as depicted in Fig. 6, for elastoviscoplastic fluids, n is less than one which renders the critical pressure gradient (or the inverse of the critical yield number; see [6,7]) to be smaller in elastoviscoplastic fluids compared to the viscoplastic fluids.
As previously observed in numerous number of studies on viscoelastic liquids, high elasticity of the fluid (i.e.high W i) results in unsteadiness, chaos and finally elastic turbulence in the inertialess regime [18]; see [55] in the context of porous media flows.In elastoviscoplastic fluid flows in a periodic cell (i.e.model porous media), we previously observed a transition to time-dependent oscillations at Bn = 10 and W i ≈ 0.5 [16], using a volume-penalization IBM approach with > 700 grid points over each cylinder surface.In elastoviscoplastic fluids, the elastic behaviour is characterised by W i × Bn as evidenced both computationally and experimentally in our previous studies [12,17,21,22].Hence, even at relatively small Weissenberg numbers, if the Bingham number is sufficiently high, the unsteadiness appears; here we quantified this behaviour in Figs.7 and 8 where the pressure drop is shown as a function of time.Also, the effect of porosity has been discussed: at low solid volume fractions (e.g.ϕ = 10%), the flow is steady at low Bn and W i, while increasing the volume fraction to ϕ = 20% or 30% makes the flow unsteady at the same low W i and Bn.In summary, the pressure drop is smaller and time-independent (steady) when the medium is more porous and the flow is at low Bn and W i numbers.Conversely, for larger values of the Bingham number, Weissenberg number and the solid volume fraction, the pressure drop becomes considerably time-dependent.This is in line with studies of viscoelastic fluids in porous media, see e.g.[56].Recent studies [57] have suggested that the number of stagnation points in the flow affects the onset of elastic turbulence.The number of stagnation points naturally increases with the solid volume fraction, but also depends on the precise geometric configuration in the disordered porous medium [57,58].Elastic turbulence results in a sudden increase in drag above a certain Weissenberg number threshold in viscoelastic flows through porous media [59].Hypothetically, if this were to happen in elastoviscoplastic flows as well when W i × Bn becomes very large, then our conclusion of drag reduction by elasticity may change.Characterising the elastoviscoplastic turbulence in porous media is, however, beyond the scope of the present work.
Due to the high cost of numerical simulations of elastoviscoplastic fluid flows, here we focused on unravelling general features of this kind of fluids in porous media flows.Future studies are required to delve into the details and also the bulk transport properties which hopefully can yield general Darcy law for such kind of fluids.Thus, we only scratched the surface here which can be a baseline for future studies to extend our work in various directions, for example, focusing on more realistic topologies (e.g.digital rocks), more sophisticated rheological models (e.g.thixo-elasto-viscoplastic behaviours), or inertial regimes.Also, in industrial applications, the porosity of the medium can be very small such as "digital rocks" but for 2D randomised monodispersed obstacles, we cannot reach that limit.This remains for future studies to be addressed.

Figure
Figure Pressure gradient versus the Bingham number: solid lines & circle symbols represent the viscoplastic simulations; dashed lines & asterisks the cases of W i = 0.05; dashed dotted lines & diamonds the cases of W i = 0.5; dotted line & squares the case of W i = 10.The black colour stands for ϕ = 10%, the purple for ϕ = 20%, and the red for ϕ = 30%.The linear scaling is marked by green dashed line for reference.

Figure 8 :
Figure 8: Instantaneous pressure gradient (∆P /L versus time) for the extreme case W i = 10 & Bn = 500 for different solid volume fractions.

Table 1 :
Drag coefficient of a cylinder in an infinite unidirectional viscoelastic fluid flow.Please note that W i and Re in this table are defined based on the particle diameter, whereas the previous studies defined it based on the radius.