A closure model for the drag coefficient of a sphere translating in a viscoelastic fluid

In many large-scale industrial applications dealing with particle-laden viscoelastic fluids, the ensemble-averaged behavior of the mixture is of most interest. The first step to parametrize this behavior is to develop an accurate expression to rapidly evaluate the drag coefficient over a broad range of kinematic parameters. The drag coefficient of a spherical particle translating in a viscoelastic matrix is strongly affected by the viscoelasticity of the fluid. In this study, we aim to parametrize the effects of fluid elasticity, especially the relaxation and retardation times, as well as inertia on the drag coefficient of a sphere translating in a viscoelastic fluid described by the Oldroyd-B model. To this end, we employed three-dimensional direct numerical simulations of viscoelastic flow past a stationary sphere. The accuracy of the numerical formulation is thoroughly tested against a number of benchmark problems consisting of steady flow past a sphere in a bounded circular or square domain filled with either a Newtonian or viscoelastic fluid. Initially, the numerical computations for the drag coefficient over a wide range of geometric and flow parameters are validated by comparison with existing data and drag correction models from the literature. The drag coefficient correction is then evaluated for unconfined flow past a sphere at different Reynolds number, Re, over a wide range of Deborah number, De < 9, and polymer viscosity ratio, 0 < zeta < 1. For small Deborah number (De < 1), the drag coefficient decreases with respect to the Stokes drag coefficient, whereas, at large Deborah number (De > 1), the drag is enhanced due to the large elastic stresses that develop on both the surface and wake of the sphere. These canonical behaviors, observed in the inertia-less flow regime (Re <= 1) are amplified as the polymer viscosity ratio approaches unity. At higher Reynolds numbers (Re > 1), the drag coefficient correction arising from viscoelasticity is found to be always bigger than unity, but smaller than the enhancement calculated in the creeping flow limit. In this regime, the formation of an axisymmetric recirculating eddy expands the wake region behind the sphere, and shifts the maximum elastic stresses away from the rear stagnation point towards the flow separation points. Motivated by the self-similarity observed at Re < 1, we propose an approximate model for the viscoelastic drag coefficient based on rational functions, using the lowest possible degree of polynomial functions required to fit the computational data with 95% accuracy, for De <= 5 and 0 < zeta < 1. This expression is a key step towards formulating a momentum-exchange model for the flow of dilute suspensions of solid particles in a viscoelastic matrix.