  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
$\langle F \rangle$ on the volume fraction of the dispersed solid phase and on
the rheology of the complex fluid matrix. Here we use direct numerical
simulations (DNS) to study the creeping flow ($Re\ll 1$) 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$ different systems, in which the normalized average fluid-particle force
$\langle F \rangle$ is obtained as a function of the volume fraction $\phi$ $(0
< \phi \leq 0.2)$ of the dispersed solid phase and the Weissenberg number $Wi$
$(0 \leq Wi \leq 4)$. From these predictions a closure law $\langle F
\rangle(Wi,\phi)$ for the drag force is derived for the quasi-linear Oldroyd-B
viscoelastic fluid model which is, on average, within $5.7\%$ of the DNS
results. Additionally, a flow solver able to couple Eulerian and Lagrangian
phases 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, in order to assess the accuracy and robustness of the
newly-developed approach for handling particle-laden viscoelastic flow
configurations with $O(10^5-10^6)$ rigid spheres that are representative of
hydraulic fracture operations.
