Finite volume simulations of particle-laden viscoelastic fluid flows: application to hydraulic fracture processes

Accurately resolving the coupled momentum transfer between the liquid and solid phases of complex fluids is a fundamental problem in multiphase transport processes, such as hydraulic fracture operations. Specifically we need to characterize the dependence of the normalized average fluid–particle force ⟨F⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle F \rangle$$\end{document} on the volume fraction of the dispersed solid phase and on the rheology of the complex fluid matrix, parameterized through the Weissenberg number Wi measuring the relative magnitude of elastic to viscous stresses in the fluid. Here we use direct numerical simulations (DNS) to study the creeping flow (Re≪1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Re\ll 1$$\end{document}) of viscoelastic fluids through static random arrays of monodisperse spherical particles using a finite volume Navier–Stokes/Cauchy momentum solver. The numerical study consists of N=150\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=150$$\end{document} different systems, in which the normalized average fluid–particle force ⟨F⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle F \rangle$$\end{document} is obtained as a function of the volume fraction ϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi$$\end{document}(0<ϕ≤0.2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(0 < \phi \le 0.2)$$\end{document} of the dispersed solid phase and the Weissenberg number Wi(0≤Wi≤4)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(0 \le Wi \le 4)$$\end{document}. From these predictions a closure law ⟨F(ϕ,Wi)⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle F(\phi,Wi) \rangle$$\end{document} for the drag force is derived for the quasi-linear Oldroyd-B viscoelastic fluid model (with fixed retardation ratio β=0.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta = 0.5$$\end{document}) which is, on average, within 5.7%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5.7\%$$\end{document} of the DNS results. In addition, a flow solver able to couple Eulerian and Lagrangian phases (in which the particulate phase is modeled by the discrete particle method (DPM)) is developed, which incorporates the viscoelastic nature of the continuum phase and the closed-form drag law. Two case studies were simulated using this solver, to assess the accuracy and robustness of the newly developed approach for handling particle-laden viscoelastic flow configurations with O(105-106)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(10^5-10^6)$$\end{document} rigid spheres that are representative of hydraulic fracture operations. Three-dimensional settling processes in a Newtonian fluid and in a quasi-linear Oldroyd-B viscoelastic fluid are both investigated using a rectangular channel and an annular pipe domain. Good agreement is obtained for the particle distribution measured in a Newtonian fluid, when comparing numerical results with experimental data. For the cases in which the continuous fluid phase is viscoelastic we compute the evolution in the velocity fields and predicted particle distributions are presented at different elasticity numbers 0≤El≤30\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\le El \le 30$$\end{document} (where El=Wi/Re\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$El=Wi/Re$$\end{document}) and for different suspension particle volume fractions.


Introduction
Understanding the force balance that governs the migration of rigid particles suspended in a viscoelastic fluid is fundamental to a wide range of engineering and technology applications.
Examples include polymer processing of highly-filled viscoelastic melts and elastomers [1], the processing of semi-solid conductive flow battery slurries [2], the flow-induced migration of circulating cancer cells in biopolymeric media such as blood [3], magma eruption dynamics [4], and hydraulic fracturing operations using solids-filled muds, slurries and foams [5,6].
In many ways, developing robust and accurate tools to simulate such behaviors may be viewed as an unsolved grand challenge in the dynamics of complex fluids, involving the effects of nonlinear material rheology, fluid inertia, elasticity, flow-unsteadiness plus many-body interactions. In a fluid with Newtonian or non-Newtonian rheology, the presence of a cloud of particles dramatically changes the transmission of stress between both phases (fluid and particles), specifically in terms of the rate at which the constituents of the mixture exchange momentum, known as hindrance effect [7,8] . When particles are suspended in a viscoelastic fluid (e.g. a polymer solution or polymer melt), the problem becomes even more challenging, because the fluid may shear-thin or shear-thicken as well as exhibit viscoelasticity and yield stress attributes [9,10]. Therefore, understanding and predicting both the bulk/macro-scale and particle-level response of these complex multiphase suspensions remains an open problem. As a first step, quantifying the momentum exchange between the constituent phases (i.e. a viscoelastic fluid matrix and a suspended phase of rigid spherical particles) remains a challenging and important problem to be solved in non-Newtonian fluid dynamics.
Because of the linear response between stress and deformation rate, the hydrodynamic behavior of rigid spheres in a Newtonian fluid has received considerable attention since the pioneering work of Stokes [11] (see for example the monographs by Happel and Brenner [12], Kim and Karrila [13] and Guazzelli and Morris [14]). In the limit of infinite dilution, and when inertial effects can be neglected, the drag force, F d , exerted by the fluid on the solid object takes the Stokes-Einstein form F d = 6πaη 0 u, where a is the radius of the spherical particle, η 0 is the fluid shear viscosity and u is the superficial fluid velocity, defined as the fluid velocity averaged over the total volume of the system [15]. Higher order corrections to the Stokes-Einstein drag acting on a single particle arising from the presence of neighboring particles have been evaluated in terms of an expansion in the particle volume fraction φ. For small packing fractions (φ < 0.10), the first few terms can be worked out analytically [16], but for larger packing fractions, the drag force has to be estimated from approximate theoretical methods [17,18], or from empirical data via experimental measurements [19]. Additionally, numerical simulations [15,20,21] can provide data to derive drag force expressions for the creeping flow of random arrays of spheres surrounded by a Newtonian fluid.
For the creeping flow of a sphere through an unbounded viscoelastic fluid, measurements of the changes to the force acting on the sphere are typically represented in terms of a dimensionless drag correction factor X(W i), which is the ratio of the measured drag coefficient compared to the well-known Stokes drag, X(W i) ≡ C D (W i)/C D (W i = 0) = C D (W i)/(24/Re) where Re and W i are the dimensionless Reynolds and Weissenberg numbers, respectively. For the inertialess flow of a viscoelastic fluid, perturbation solutions predict the departure of the drag from the Stokes result to be quadratic in the Weissenberg number for spheres [22], and linear in the case of long rod-like particles [23]. Recently, two reviews comparing experimental data with computations for non-Brownian suspension rheology with non-Newtonian matrices have been published [9,10]. Tanner [10] compared and contrasted inelastic fluids with rate-dependent viscosity, materials with a yield stress, as well as viscoelastic fluids -highlighting the need for rheological modeling improvement, possibly with multiple relaxation times. Additionally, he concludes that several aspects of suspension rheology, such as roughness, ionic forces, particle shape, and polydispersity all need to be addressed. Finally, Tanner [10] also reported experimental results for steady viscometric flows, unsteady shear flows and uniaxial elongational flows. However, good agreement between computation and experiment is scarce, because there are, as yet, few computational studies which allow careful comparison with experimental data, further emphasizing that progress in rheological modelling and improved computational methods are needed. In a recent perspective, Shaqfeh [9] notes that the foundations for the development of suspension mechanics in viscoelastic fluids, as well as the development of computational methods to accurately simulate with particle-level non-Brownian suspensions, have been established. Nevertheless, numerous unanswered questions remain, including the rheological behavior of these suspensions for different matrix fluid rheologies, particle shapes, deformability, flow histories, etc. All these questions can be addressed, in principle, by employing theoretical/computational frameworks to systematically explore the coupling between the kinematics and momentum distribution of the fluid phase and the resulting evolution of the dispersed particulate phase. Jain and Shaqfeh [24] performed 3D transient simulations of the bulk shear rheology of particle suspensions in Boger fluids for a range of W i ≤ 6 and finite strains and calculated the per-particle extra viscosity of the suspension. They categorize the per-particle viscosity calculations as contributions from either the particle-induced fluid stress (PIFS) or stresslet contributions. It was concluded that in the dilute limit, the PIFS increases monotonically with shear strain; however, the stresslet contribution shows a non-monotonic evolution to steady state at large W i. The total combined per-particle viscosity contribution, however, shows a monotonic evolution to steady state. Additionally, Jain and Shaqfeh [24] performed multiple-particle simulations using the IB method to examine the effect of particleparticle hydrodynamic interactions on the per-particle viscosity calculation. It was concluded from transient immersed boundary simulations that the steady values of per-particle viscosity increase with φ, but the per-particle contribution to the primary normal coefficient was independent of φ (up to 10% particle volume fraction) at the two values of Weissenberg number investigated (W i = 3 and 6).
The nonlinear interactions of fluid inertia with viscosity and elasticity cause unexpected phenomena (e.g. negative wakes, shear-induced migration/chaining) in the dynamics of particles suspended in a viscoelastic matrix [25,26,27]. These interactions may be expected to change the evolution in the viscoelastic drag correction factor with Weissenberg number.
Extensive research efforts over the past 20 years have been directed at the elucidation of the role of fluid rheology and wall effects on the drag of a sphere and the wake developed behind it. Excellent reviews are available in the literature [28,29,30]. There have been a number of computational studies investigating the effect of fluid rheology on the motion of the sphere; however, a common limitation is that results are not convergent for Weissenberg numbers beyond a certain critical limiting value (i.e. typically W i c ≈ 2 or 3) [28]. In the work we extend viscoelastic suspension flow calculations up to W i = 4, by employing the log-conformation approach [31,32,33,34].
Our previous work [35] proposed, for the first time, a closure model for the viscoelastic drag coefficient of a single sphere translating in a quasi-linear viscoelastic fluid, which can be well described by the Oldroyd-B constitutive equation. In the present work, we extend the proposed model in order to be able to describe moderate volume fraction viscoelastic suspensions (φ ≤ 0.2), which are commonly encountered in a wide range of industrial operations. We focus on non-colloidal suspensions with Newtonian and viscoelastic fluid matrices, and the net effect of other particles in the flow is studied by computing an effective average drag force acting on a particle [36].
For Newtonian matrices, Brinkman [17] presented a modification of Darcy's equation in porous media, in which the viscous force exerted on a dense suspension of rigid particles by the Newtonian fluid is calculated. The main idea is that, from the point of view of a single particle, the other distributed particles effectively act as a porous or Darcy medium. Durlofsky and Brady [37] performed numerical simulations to compare against the predictions of Brinkman's model; however, they only obtained good agreement with the analytical solution of Brinkman for very small volume fractions up to φ = 5%, emphasizing the importance of including the configurational effects of more distant particles. The reason that the Brinkman approach starts to break down at what appears to be a rather low volume fraction is related to the fact that at φ = 5% a characteristic inter-particle spacing is only slightly larger than four particle radii; so that neighboring particles are in fact quite close together and hydrodynamically interact. For dense random arrays of spheres (φ ≥ 40%) the empirical Carman-Kozeny (CK) relation [19] is found to describe well the drag force exerted by the fluid flow in the dispersed phase. The idea behind the CK relation is that the suspended medium can be considered as a system of tortuous channels, from which the pressure drop across the porous medium is calculated using the Darcy equation [38]. In the present work, we report results for volume fractions of φ = 4%, 8%, 12%, 16% and 20%, representative of semi-dilute non-colloidal suspension behaviour [39].
Additionally, we note that fully-resolved particle-laden viscoelastic solvers [40,41] are presently only able to directly resolve O(10 3 ) particles [42], which limits their application to large industrial case studies [43]. To overcome this limitation we implement an Eulerian-Lagrangian viscoelastic solver (DP M viscoelastic), which employs the closure drag model that we develop here for moderately dense suspensions with a viscoelastic matrix fluid to quantify the momentum exchange between the two constituent phases (a moderate volume fraction of rigid spherical particles and a non-shear thinning viscoelastic matrix fluid).
The present paper describes the simulation method employed to measure the drag force on randomly-dispersed particle arrays immersed in viscoelastic fluids that can be described by the quasi-linear Oldroyd-B constitutive equation, which predicts constant values of the shear viscosity and first normal stress coefficient. The extension of this work to more complex viscoelastic matrix-based fluids, for example models predicting shear-thinning fluid behavior (such as the Giesekus model), is also currently being studied. In this case the magnitude of the stresses is typically smaller and easier to resolve computationally, but the dimensionality of the problem is higher due to the additional nonlinear model parameter(s) required to characterize the shear-thinning. To accomplish the required developments, an open-source library, OpenFOAM [44], was modified to be able to calculate the average drag force acting on random particle arrays. The latter information is then used to formulate a new drag force correlation for the creeping flow of an Oldroyd-B fluid through randomly distributed arrays of spherical particles, with solid volume fractions 0 < φ ≤ 0.2, over a range of Weissenberg numbers (W i ≤ 4). To ensure stability at high W i, the polymer stress contribution is computed using the log-conformation formulation [31,32,33,34]. Finally, to the best of the authors' knowledge, the current work presents for the first time an Eulerian-Lagrangian solver, in which the fluid continuum phase has viscoelastic rheology and the dynamics of the particulate phase are computed following a discrete particle method (DPM). The momentum transfer between both phases is calculated using the drag force law proposed in the present work (see Section 4.2).
The newly-developed DP M viscoelastic solver is employed to predict particle settling effects in rectangular channels (which mimics a vertical fissure such as those encountered in gas/oil extractions during hydraulic fracturing operation) and in an annular pipe (a model of pumped transport of a suspension along a drill string during horizontal drilling operations).
The paper is organized in the following manner: in Section 2 we present the governing equations and numerical method used to compute the solution of the appropriate balance equations for the viscoelastic fluid flows considered in this work. Section 3 provides the details of the simulation methodology used to compute the drag force values exerted by the fluid on the particulate phase. In Section 4 these results are verified for Stokes flow of a particle suspension dispersed in a viscous Newtonian fluid. The results are then used to derive a new closure law for the average fluid-particle drag force acting on random arrays of spheres immersed in an Oldroyd-B fluid. Section 5 is dedicated to description of the development of the DP M viscoelastic solver for up-scaled three-dimensional simulations of particle-laden viscoelastic flows, in which the dispersed phase is modeled by the discrete particle method.
Additionally, we illustrate the capability of the newly-developed code to solve challenging physical problems, specifically two canonical proppant transport problems, which are commonly encountered in hydraulic fracturing operations. Finally, in Section 6, we summarize the main conclusions of this work.

Governing equations
Following the work of Faroughi et al. [35], as a first step, we consider the problem of moderately dense (0 < φ ≤ 0.2) suspensions constituted from a continuum viscoelastic matrix fluid and a monodisperse static random array of rigid spheres. For the continuum fluid phase the familiar Oldroyd-B constitutive model was chosen, representing an elastic fluid with a constant shear viscosity, which has been shown by Dai and Tanner [45] to fairly well describe the response of highly elastic Boger fluid suspensions in steady shear and uniaxial elongation. By adopting the Oldroyd-B model, we confine the dimensionality of the viscoelastic fluid calculations to only two degrees of freedom (i.e. the relaxation and retardation times, or equivalently).
However, the consideration of moderately dense suspensions increases the dimensionality of the problem to four degrees of freedom, due to the addition of two more variables; the particle volume fraction present in the suspension and the number of random particle configurations studied (to obtain statistical significance in the DNS results). Thus, for the current study, we have also fixed the retardation time of the fluid at β = η S /(η S + η P ) = η S /η 0 = 1/2, where η 0 is the total matrix fluid viscosity, with η S and η P being the solvent and polymeric viscosities, respectively. Within these constraints, the drag correction expression developed in this study will form a foundation for higher-dimensional parameterizations, which should also consider a range of solvent viscosities as well the effect of more complex fluid rheology (e.g. shear-thinning) on random arrays of particles in viscoelastic fluids, by using machine learning algorithms such as convolutional neural networks to capture the non-linear effects of all constitutive parameters on the resulting drag coefficient expressions acting on the particle arrays.
The dimensionless conservation equations governing transient, incompressible and isothermal laminar flow of an Oldroyd-B fluid are given by where the following dimensionless quantities are used with L and U being the characteristic length and velocity values, respectively, x the position vector, u the velocity vector, t the time, p the pressure and τ P the polymeric contribution to the extra-stress tensor. As mentioned above, the retardation ratio is fixed with the value of β = 0.5 for all the calculations performed in this work.
For the present problem with L = a, where a is the radius of a single suspended particle, and U the average fluid velocity at the inlet of the channel, we define the Reynolds and Weissenberg numbers as follows, where ρ is the fluid density and λ is the relaxation time. Notice that for the case of a Newtonian fluid flow λ = 0 and η 0 = η S .
To ensure computational stability over a wide range of fluid elasticities, including suspension flows at high Weissenberg number, we incorporate the log-conformation approach for calculating the polymeric extra-stress tensor. In the present work, we follow the implementation of the log-conformation approach in the OpenFOAM computational library [44], presented in Habla et al. [33] and Pimenta and Alves [34]. Details on the mathematical for-mulation behind the log-conformation approach can be found in the original works of Fattal and Kupferman [31,32].

Numerical method
The equations presented in Section 2.1 are discretized using the finite-volume method (FVM) implemented in the OpenFOAM framework [44].
Pressure-velocity coupling was accomplished using segregated methods, in which the continuity equation is used to formulate an equation for the pressure, using a semi-discretized form of Eq. (1) [46]. The resulting equation set is solved by a segregated approach, using the SIMPLEC (Semi-Implicit Method for Pressure-Linked Equations-Consistent) algorithm [47], which does not require under-relaxation of pressure and velocity (except for non-orthogonal grids, where the pressure needs to be under-relaxed [34]). Additionally, the computational cost per iteration of this algorithm is lower than in the PISO (Pressure-Implicit Split Operator) algorithm [48], because the pressure equation is only solved once per cycle. The coupling between stress and velocity fields is established using a special second-order derivative of the velocity field in the explicit diffusive term added by the iBSD (improved both-sides diffusion) technique [49]. The velocity gradient is calculated using a second-order accurate least-squares approach, and the diffusive term in the momentum balance is discretized using second-order accurate linear interpolation. For non-orthogonal meshes the minimum correction approach is used, as explained in Jasak [50], in order to retain second-order accuracy. The advective terms in the momentum and constitutive equations are discretized using the high-resolution scheme CUBISTA [51] following a component-wise and deferred correction approach, enhancing the numerical stability. The time derivatives are discretized with the bounded second-order implicit Crank-Nicolson scheme [52]. Here, a Poisson-type equation for the pressure field is solved with a conjugate gradient method with a Cholesky preconditioner, and the linear systems of equations for the velocity and stress are solved using BiCGstab with Incomplete Lower-Upper (ILU) preconditioning [53,54,55]. The absolute tolerance for pressure, velocity and stress fields was set as 10 −10 . The simulations are performed including transient terms, but the time marching is used only for relaxation purposes as we will just be looking for the steady-state solution, i.e., when the drag coefficient ceased to vary in the third decimal place.

Simulation methodology
In order to develop our computational methodology, we address only non-colloidal suspensions with viscoelastic matrices, and focus on both dilute and moderately dense suspensions.
Following Housiadas and Tanner [36], we consider the effect of other particles in the flow by assuming that from the point of view of a single particle at any instant, the remaining particles act effectively like a porous medium [17].
is as close as possible to the desired packing fraction, where n s is the number of spheres enclosed in the square duct (whose volume is LH 2 ) and a is the sphere radius. For the purpose of simulating the proppant transport phenomena that occur during hydraulic fracturing operations, we consider in this work moderately dense suspensions where the particle volume fraction range is 0 < φ ≤ 0.2 [5,43,56]. Following Hill et al. [20] we note that the system to be studied must include a sufficient number of spheres, n s , to minimize artifacts and statistical oscillations coming from the finite size of the computational domain.
In practice, n s is chosen to be large enough to avoid periodic artifacts (typically 24 ≤ n s ≤ 122, see Table 1 in Section 4.1), and statistical uncertainty is reduced by ensemble averaging the results from n c random sphere configurations (in this work n c = 5 was found to be sufficiently large to obtain a standard error for the average below 5% of its actual value, as shown in Section 4, Tables 1 and 2).
The numerical model employed in this work was comprehensively tested against a similar computational challenge (see Faroughi et al. [35] for the bounded and unbounded flow past a single sphere in a fluid described by the Oldroyd-B constitutive model). The meshes employed in this work have the same level of mesh refinement as the most refined mesh (M1) used by Faroughi et al. [35], which resulted from a grid refinement study.
In all cases, the magnitude of the fluid velocity imposed at the inflow was such that the Reynolds number based on the particle diameter D, Eq. (5a), was equal to Re D = 0.05, representative of creeping flow conditions. At this point we should note that there exists some ambiguity in the literature [15] on the proper definition of the drag force, specifically, if the pressure gradient term should contribute to the drag force or not. It is known that the two definitions differ by a factor (1 − φ), i.e., the relation between the total average force that the fluid exerts on each particle, F t , and the drag force, F d , which results from the friction between the particle and the fluid at the surface of the particle, is in some literature (Hill et al. [20]), the total force on the particle is defined as the drag force.
In this work, the results will be presented in terms of an average dimensionless drag force, F , which is defined as where u = (1 − φ)U is the superficial fluid velocity, F t,i is the average drag force on the random array of spheres computed as F t,i = 1 ns ns i=0 F t,i , where F t,i is the drag force on sphere i from an ensemble of n s spheres, and e x is the unit vector in the x-direction. The denominator on the right-hand side of Eq. (6) is the Stokes drag force, obtained in the limit of infinite dilution and when inertial effects can be neglected [11]. In this work, the uncertainty in the computed average force, referred as the standard error, is calculated from Note that the factor n c − 1 in the denominator on the right-hand side of Eq. (7) corrects for the fact that there are n c − 1 degrees of freedom, since the average is used to calculate the variance in the numerator [57]. This is important since the number of random configurations, n c , used to calculate the average of F is small.

Results
: simulation of flow in random particle arrays

Verification: Stokes flow of suspensions with Newtonian fluid matrices
In this section the creeping flow of random arrays of spheres surrounded by a Newtonian fluid is studied. This way the simulation methodology presented in Section 3 can be verified against results found in the literature.
One of the earliest drag force models for describing Stokes flow through an array of spherical particles is the Carman [19] relation, This relation is only valid for dense arrays ((1 − φ) 1), which can be seen by the fact that it does not have the correct limit F → 1 for φ → 0. For the limit of dilute systems, Kim and Russel [18] derived a closed-form expression for F : The computational results obtained by Hill et al. [20], using Lattice-Boltzmann simulations, were found to be in very good agreement with Kim and Russel's [18] drag force expression (Eq. (9)) for dilute arrays of particles (φ ≤ 0.1).
Subsequently, several expressions have been developed to find an accurate drag force model that is valid over the full solid fraction range. Using a modification of the Darcy equation [38], Brinkman [17] derived the well-known drag force model, Koch and Sangani [21] proposed the following expression for the drag force: which for low solid volume fraction is equal to Eq. (9) to O(φ ln φ), whereas for large solid volume fractions is the drag force given by the Carman expression Eq. (8). Finally, van der Hoef et al. [15] presented a best fit to simulation data obtained using a Lattice-Boltzmann method (for φ ≤ 0.6), which takes the following simple form: which is the Carman expression with a correction term for the limiting case of φ → 0. Recently, Faroughi and Huber [7] predicted the correction to the drag coefficient for monosized spherical particles as: where φ m denotes the maximum random close packing fraction, which, for monosized spherical particles, we take to be φ m ≈ 0.637 [58], and β is a geometrical proportionality constant, which is related to the shape of the streamtube in the real flow field. The value of β = 0.65 provides the best fit to numerical simulations shown in Fig. 2 (8)), Kim and Russel [18] (Eq. (9)), Brinkman [17] (Eq. (10)), Koch and Sangani [21] (Eq. (11)), van der Hoef et al. [15] (Eq. (12)) and Faroughi and Huber [7] (Eq. (13)) expressions are also represented in Fig. 2 for comparison. For all the simulated particle volume fractions, the dimensionless drag force F obtained from our numerical algorithm is similar  Table 1 for φ = 0.2). For a more detailed discussion regarding the excluded volume provided by impenetrable channel side-walls refer to Appendix A. The sphere volume fractions, the number of spheres used on each simulation, the number of mesh points, the dimensionless average drag force F on the spheres and the respective standard errors are listed in Table 1. In all cases, the standard errors of F , which measure the statistical accuracy achieved by averaging the results with five random configurations, were below 3.5% of the average. In Fig. 3 we show normalized axial velocity contours obtained from the numerical simulation of random particle arrays in a channel filled with a Newtonian fluid. As can be seen from the velocity contours, as we increase the particle volume fraction more of the fluid is forced to flow through the tortuous paths in the interstitial spaces between the spheres, rather than as a continuous fluid stream that is mildly perturbed by widely separated spheres. This is also visible by the higher magnitude of the fluid velocities near the channel walls where the fluid is squeezed. Notice that for the higher particle volume fraction employed φ = 0.2 we have also conducted simulations where particles can be located in the excluded volume region provided by the rigid bounding walls (φ * = 0.2), and the dimensionless average drag force F remains similar (approximately 2% higher) to the case with impenetrable side walls where the particles are located only in the central portion of the channel beyond an excluded volume region of thickness a that is adjacent to the bounding side walls region (see Fig. 2 and Table 1). Placing spheres uniformly throughout the entire domain results in a more uniform velocity profile across the channel cross-section. The dimensionless average drag force F (φ, W i) on the spheres (see Eq. (6)) and the respective standard errors are listed in Table 2 for the different kinematic conditions stated above. In all cases, the standard errors of F (φ, W i) , which measure the statistical accuracy achieved by averaging the results with five random configurations, were below 4.7% of the average drag force.   Additionally, we show in Fig. 4 The numerical results presented in Figure 5(b) show that this rescaled drag force can be considered independent from any statistically-significant trend with Weissenberg number for the range of solid volume fractions computed in this work; i.e, when we normalize the dimensionless average drag force we compute in our simulations by F 0 (W i) it helps to collapse F (φ, W i) at any given value of φ ≤ 0.20.
In the last section we obtained a good agreement between the numerical results for the average drag force exerted on an ensemble of particles immersed on a Newtonian fluid with the model given by van der Hoef et al. [15], we thus propose fitting an equation of the same form to the computational results obtained in a viscoelastic fluid, i.e., where F 0 (W i) is the infinitely dilute (φ → 0) result for the drag force presented in Faroughi et al. [35] for ζ = 0.5 (cf. Eq. (14)). Notice that we could also have used a similar expression for the correction to the drag coefficient as that given by Faroughi and Huber [7] in Eq. (13) to fit our viscoelastic results, however, due to its additional simplicity we have here chosen the functional form of the van der Hoef et al.  In Fig. 6 and Fig. 7 we show contours of the dimensionless first normal stress difference, defined as (τ xx − τ yy )/(η P U/a), obtained from the numerical simulations. The first normal stress difference is mainly generated near the no-slip surfaces and in the wake of each of the spheres. Notice that increasing the fluid elasticity (increasing W i number), i.e. from the left to the right panels in Fig. 6, promotes a strong elastic wake and an increase in the magnitude of the first normal stress difference. Only by use of the log-conformation approach for computing the polymeric extra-stress tensor components were we able to stabilize the numerical algorithm. Additionally, increasing the particle volume fraction, i.e., from top to bottom panels, increases the magnitude of the first normal stress difference that is generated on the front stagnation point of the particles. Finally, from Fig. 7 we see that the magnitude of the first normal stress difference near the rear stagnation point is much larger than the one developed upstream near the front stagnation point, and this difference becomes progressively larger for higher W i.  Fig. 7(b)).

Proppant transport during the hydraulic-fracture process
In this section we develop a computational framework, based on the Eulerian-Lagrangian formulation [60], capable of numerically describing, as a proof-of-concept, proppant transport in a viscoelastic matrix-based fluid that can be characterized by the Oldroyd-B constitutive model. The newly-developed algorithm takes into account the effect of the particle volume fraction (the Lagrangian phase) on the viscoelastic fluid phase (the Eulerian phase). For this purpose, we extend the formulation presented in Fernandes et al. [60] (and references therein) to be able to take into account the viscoelastic behavior of the fluid.
Consider the motion of an incompressible viscoelastic fluid phase in the presence of a secondary particulate phase, which is governed by the volume-averaged continuity equation and Cauchy momentum equation where f is the fluid porosity field satisfying f = 1 − φ, U f is the fluid velocity, P is the modified pressure (p/ρ f , with p being the dynamic pressure and ρ f the fluid density), g is the gravity acceleration vector and the fluid-phase stress tensor τ f is given by: where τ P is the polymeric extra-stress tensor computed using the Oldroyd-B viscoelastic matrix-based constitutive model given by Eq. (3).
The two-way coupling between the fluid phase and particles is enforced via the source term S p in the momentum balance equations, Eq. (17), of the fluid-phase. Because the fluid drag force, F d,i , acting on each particle i is known (see Section 4.2), then according to Newton's third law of motion the source term is computed as a volumetric fluid-particle interaction force given by: where V cell is the volume of a computational cell, and N p is the number of particles located in that cell.
In this work, we consider two different formulations to describe the contact between particles, the spring-dashpot model [61] and a Multi-Phase Particle In Cell (MPPIC) model [62].
The former allows us to explicitly handle each contact between two particles and, therefore, is very computationally intensive. The latter can be used to represent the particle collisions on average without resolving particle-particle interactions individually. In the MPPIC method the particle-particle interactions are computed by models which utilize mean values calculated on the Eulerian mesh [63]. For that purpose, in the present work we have employed a collision damping model to represent the mean loss in kinetic energy which occurs as particles collide, and helps to produce physically realistic scattering behaviour [63]. Finally, a collision isotropy model is also employed to spread the particles uniformly across cells [63].
Two case studies were performed to validate the newly-developed DP M viscoelastic solver.
In the first case, we study proppant transport and sedimentation in a long conduit of rectangular cross section, a typical geometry to study flow in hydraulic fracturing [43]. In the second case, we study the segregation phenomena which occurs in cement casing for horizontal wells.
Both case studies were performed using Newtonian and viscoelastic carrier fluids. For the first case, particle collisions are modelled using the MPPIC model in order to handle O(10 6 ) particles, and for the latter, the Hertzian spring-dashpot model [64] is employed with a total of 125, 000 particles. The following sections present comparisons between the numerical results obtained for the aforementioned case studies and results found in the scientific literature.

Rectangular channel flow
Despite many advances in hydrocarbon reservoir modeling and technologies [65,66,67,68,69] especially for unconventional resource development, the efficiency of hydrocarbon recovery in shale reservoirs is still very low [70]. One of the leading issues causing this inefficiency is the lack of proper proppant placement in the fracture networks. Proppant emplacement within fractures directly impacts productivity because it controls both short-and long-term conductivities of the fractured wells [71]. Proppant particles must be carried over large distances to ensure successful placement, which requires a spatially homogeneous distribution of particles [6]. However, flows of non-Brownian particles, such as proppant, often result in non-homogeneous patterns, in which particle sedimentation is commonly observed [43]. The   Table 3 gives the pressure drop results corresponding to each of the mesh refinement levels employed at the prescribed inlet flow rate Q. Additionally, Table 3    As noted by Meeker et al. [43], the onset of nonlinear behavior is most likely to be caused by some deformation of the poly dimethylsiloxane elastomer at higher pressures. Therefore, the flow rate of our numerical tests is set at Q = 100 cm 3 /h in the following studies with suspensions. This results in an average fluid velocity U = Q/HW = 5.6 × 10 −3 m/s, which corresponds to Re W = 0.06832, confirming that we are in the creeping flow regime. We also note that U is much greater than the average particle sedimentation velocity U Stokes = 3 × 10 −5 m/s, and thus, our simulations are performed under favorable transport conditions with minimal settling at the entrance. In fact, the slope of the trajectory of a sedimenting particle being transported at this average speed, U Stokes /U ≈ 5×10 −3 , is similar to the ratio of channel height to length H/L, meaning that most of the initially suspended particles entering the channel should settle as they approach the channel exit.  This is particularly noticeable in the second and third channel observation sections, in which the sediment height increases markedly. Then, the sediment height ceases to grow further quite abruptly, although the particle suspension continues to flow through the channel. This steady-state behavior persists for the duration of the flow and represents a balance between sedimentation and shear-driven resuspension. Upon cessation of flow we immediately observe a collapse in the dense phase height, i.e. the particle phase settles further to form a more compact final sedimented state. Figure 11 shows the contours of the fluid porosity distribution . This value changes significantly based on the shape of particles [74] and the size ratio between particles in the pack [75,8]. Contours of f 0.4 thus correspond to effectively solid deposits of particles. Again we can conclude that the steady-state sediment height during suspension flow increases with the initial suspension volume fraction, and that following the cessation of flow the sediment bed compactifies and reduces in height.  To quantitatively define the steady-state sediment height, h, under flow and the static sediment height, h 0 , once the flow of the suspension has ceased, we compute the average fluid porosity¯ f over a local section as¯ and define appropriate characteristic values for¯ f to quantify h and h 0 . For particular flow conditions of a non-Brownian suspension flowing at Q = 100 cm 3 /h, Meeker et al. [43] observed the buildup of a dense but flowing sediment that rapidly reaches a steady-state height h. The existence of this steady-state flowing sediment implies that the proppant flux leaving the channel equals that entering the channel, and thus, an "efficient" proppant transport occurs. Knowing this fact, we define the criteria to compute h as¯ f = 1 − φ i (see Fig. 12(a)). Because the flow is at a low Reynolds number (Re W = 0.06832), the relevant mechanism of sediment transport must be viscous resuspension (flow of an "expanded" sediment at an equilibrium height, and its subsequent "collapse" once the flow ceases [76,77]). To quantify h 0 , we quote the work of Meeker et al. [43] stating that for quiescent conditions the packing volume fraction when water is the suspending fluid is φ p ≈ 0.58, which is close to φ m . However, when the 85:15 w/w % glycerol/water mixture (viscosity ≈ 0.1 Pa.s) was employed as the suspending fluid, the packing fraction decreased to φ p ≈ 0.5. Following this, we consider the criterion to compute h 0 as¯ c = 0.5 to represent a dense/compact suspension bed (see Fig. 12(a)).
An example of these computations is shown in Fig. 12(b), which depicts the evolution of   Fig. 13(a)) at El = 0, when analyzing the particle distribution at x = L/3 (top left images), we notice that there is a significant sedimentation layer where the particle velocity is zero. In the middle and top zones of the channel, both the matrix fluid and particles flow smoothly and axially along the channel. The particles continue to slowly sediment and eventually join the deposited layer with (U x ) p → 0. For the quasi-linear Oldroyd-B viscoelastic fluid ( Fig. 13(b)) at El = 30, when analyzing the particle distribution at x = L/3 (top right images), we notice that the distribution of particle velocities is almost uniform along the channel height, and only a thin layer of sedimentation is observed. This is in contrast to the Newtonian fluid behavior. At x/L = 2/3, the behavior of the particles and fluid constituents is not substantially changed from that at x/L = 1/3, and, therefore, there is no significant particle settling zone along the channel floor for an elastic fluid with El = 30.

Annular pipe flow
Particle segregation in pumped concrete is one of the big challenges encountered when creating casing for horizontal drilled wells [78,79]. In this case, particulate solids tend to segregate axially across the pipe due to differences in the size, density, shape and other properties of the constituent phases. The corresponding increase in the percentage of cementitious particles in the bottom part of the casing increases the chance of shrinkage and formation of cracks in the upper portion of the cemented casing. These cracks, often large in size, can easily transport hydrocarbons and other toxic chemicals into the formation which is a concern. Tuning the rheology of the conveying fluids systematically by considering the hindrance effect (i.e., reduction in the relative settling velocity of a particle due to the presence of other particles) can help minimize this issue.
Here we study numerically the particle segregation in a simplified annular pipe geometry.
The setup used to study the particle segregation is shown schematically in Fig. 15. The channel interior, which is used to mimic an horizontal well, has an annular cross section with inner The initial setup for this computational study is shown schematically in Fig. 16(a). A total of 125,000 particles, representing 1% of the annular cavity volume, is used in this case study. The particles are distributed evenly throughout the stagnant fluid at time zero (see Fig. 16(a)). The particle positions at time t = 0 are generated using a nearest neighbor algorithm. Gravity is applied vertically across the thin annular geometry (g = −ge z ), breaking the azimuthal symmetry and mimicking the onset of concrete particle settlement right after injection is stopped. The goal here is to capture the settling dynamics and test our numerical code to reproduce those dynamics. The code can be then used to analyze different rheological tuning mechanisms to minimize settling over the required time-scale for the concrete to harden.    For these two elasticity numbers the particle distributions are similar, with a settling zone and avalanche zones at the bottom and lateral walls of the annular pipe domain, respectively.
Additionally, in the settling zone, a backflow of particles occurs due to fluid displaced by the sedimentation and net accumulation of particles in this region. At the north pole (point N in Fig. 17) of the inner cylinder wall the particles have a backflow velocity, which makes them bounce and slide along the inner cylinder wall. Subsequently, the particles approach the most unsteady settling zone of the annular pipe channel, where a mixture of fluid backflow and which indicates a migration of the particles to the avalanche zone. In fact, from the particle velocity distributions, we see that the stronger migration of the particles to the avalanche zone at El = 0.1 causes an increase in the suspension bed height when comparing to the higher elastic case with El = 5. In fact, the calculated final packed bed height for the cases where El = 0 or El = 0.1 is 4.5 mm and for El = 5 is 3.5 mm.

Conclusions
Direct numerical simulations (DNS) of random arrays of spherical particles immersed in Newtonian and constant-viscosity viscoelastic fluids were performed using a finite-volume method. The overall procedure solves the equations of motion coupled with the viscoelastic Oldroyd-B constitutive equation using a log-conformation approach, with a SIMPLEC (Semi-Implicit Method for Pressure-Linked Equations-Consistent) method. The drag forces on individual particles were calculated with the aim of providing an approximate closed-form model to describe numerical simulation data obtained for the unbounded flow of Newtonian and Oldroyd-B fluid past random arrays of spheres. This expression can then be integrated into a Eulerian-Lagrangian solver that enables coupled simulations of the fluid flow and particle migration over a wide range of kinematic conditions. For this purpose, the DNS consisted of a total of 150 different configurations, in which the average fluid-particle drag force is obtained for solid volume fractions φ (0 < φ ≤ 0.2) and Weissenberg number W i (0 ≤ W i ≤ 4).
The proposed DNS methodology was first tested and verified for the creeping flow of random arrays of spheres immersed in a Newtonian fluid. It was found that the numerical results obtained agree with the Lattice-Boltzmann results of Hill et al. [20] and can be described by the best-fit model of van der Hoef et al. [15]. Statistical accuracy was achieved by averaging the DNS results at each value of φ over five random configurations, resulting in errors below 3.5% of the average drag force. Subsequently, the same DNS methodology was used to per- employed together with a viscoelastic constitutive equation to describe the fluid flow, and a discrete particle method is used to update the particle movements. This approach guarantees the coupling between the dynamics of the continuous fluid and the discrete solid phases, by imposing a two-way coupling between the two phases. The coupling is provided by momentum transfer through the drag force expression proposed here, which is exerted by the fluid on the solid particles. Additionally, we consider two different formulations to describe the contact between particles, the Hertzian spring-dashpot and Multi-Phase Particle In Cell (MPPIC) models.
As a proof-of-concept, the newly-developed algorithm was assessed for accuracy in two case studies. First, we studied the proppant transport and sedimentation during pumping (a phenomenon typical of hydraulic fracturing operations) in a long channel of rectangular cross section. For the case in which the fluid matrix is Newtonian, the resulting axial distribution of particle sedimentation profiles was compared with experimental data available in the literature for suspensions formulated with Newtonian matrix fluids and different initial particle volume fraction, and good agreement was obtained. Subsequently, the DP M viscoelastic solver was tested on the same problem using an Oldroyd-B fluid. Analysis of the particle distribution and fluid velocity profiles at an elasticity number of El = 30, showed that fluid elasticity inhibits the rate of particle settling and prevents the formation of a dense sedimented layer along the floor of the channel.
Subsequently, segregation phenomena which occurs when pumping a casing material along horizontal wells was also studied in an annular pipe domain. Numerical simulations using a Newtonian fluid were performed, and we were able to capture the avalanche and dome build-up effects observed in experimental observations of the particle distributions [79]. Additionally, a viscoelastic fluid was also employed at two different elasticity numbers El = 0.1 and 5. The particles were found to sediment with two markedly contrasting zones, a highly disordered and unsteady region where a mixture of fluid backflow and gravity-induced settling velocities are present and a sedimented zone where particles are closely packed together and the fluid velocity is almost zero. It was found that the stronger migration of the particles to the avalanche zone at El = 0.1 cause an increase in the suspension bed height when comparing to the higher elastic case with El = 5.
In summary, the DNS computational methodology presented here allows us to construct

Appendix A
In the DNS study presented in section 4, we considered two different domain configurations, one with spheres having centroids in a wall region of thickness a around all four lateral edges of the flow domain and another in which the sphere centroids are excluded from this wall region. We refer to these cases as the no-excluded volume and excluded volume configurations, respectively. As shown in Fig. 18(a) when spheres are allowed to be located in the wall region (blue color), i.e., when their centroid is located less than one radius from the wall, then the boundary acts as a perfectly periodic wall. In the opposite case the boundary walls exclude the spheres and act like rigid stress free periodic walls. In fact, based on Fig. 18(b), we can calculate the probability of a single sphere being located in the wall region. Assuming a square cross-section with a width of 8a, the total cross-sectional area is 64a 2 . Regarding the blue annular area, i.e., the region which excludes the spheres near the wall, the area is of (64a 2 − 36a 2 = 28a 2 ). Hence, the probability of a randomly placed sphere being located in the excluded region is equal to 28a 2 /64a 2 = 0.4375 and the overall/fraction area of spheres (of volume fraction φ) in this region can be as large as 0.4375φ. Therefore, as φ increases it is important that when particles are randomly distributed in the domain they should be allowed to be placed with centroids near the walls. In Fig. 19, we show contours of the velocity magnitude in the transverse y − z plane for configurations with an excluded volume and no-excluded volume region. From the distribution of velocity magnitude contours, it can be seen that in the excluded volume configuration (i.e., Fig. 19(a)), the larger local concentration of rigid impenetrable spheres in the middle of the square channel push the strongest fluid flow ownwards towards the walls causing a stagnant region near the channel center. On the other side, in the no-excluded volume configuration (i.e., Fig. 19(b)), the fluid flow is more evenly distributed across the entire channel. This affects the average drag force exerted on the spheres as shown in Table 1 and Fig. 3. : Steady flow field around one representative random array of particles in a channel filled with Newtonian fluid. Contours of the velocity magnitude field u (with the inflow direction pointed out of the plane of the page) are represented for particle volume fraction φ = 0.2, in the y − z plane with (a) excluded volume near the walls and (b) a no-excluded volume configuration. In the latter case the velocity field is more evenly distributed across the entire cross-section of the domain.