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 > 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) 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 <F> is obtained as a function of the volume fraction φ (0 < φ ≤ 0.2) of the dispersed solid phase and the Weissenberg number Wi (0 ≤ Wi ≤ 4). From these predictions a closure law < F >( Wi,φ ) for the drag force is derived for the quasi-linear Oldroyd-B viscoelastic fluid model (with fixed retardation ratio β = 0.5) which is, on average, within 5.7% of the DNS results. Additionally, 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, 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. 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 (where El = Wi/Re ) and for different suspension particle volume fractions.