Increase of turbulent drag by polymers in particle suspensions

We study the effect of spherical particles on the turbulent flow of a viscoelastic fluid and find that the drag reducing effect of polymer additives is completely lost for semidense suspensions, with the drag increasing more than for suspensions in Newtonian fluids. This different behavior is due to three separate effects. First, polymer stretching is reduced by the presence of rigid particles, thus canceling the drag reducing benefit of the viscoelastic fluid. Second, drag increase is provided by the growth of the particle and polymeric shear stresses with the particles, due to larger shear rates in the vicinity of the particle surface. Third, particles migrate towards the wall due to the shear-thinning property of the fluid, thus enhancing the particle near-wall layer and further increasing the drag.

drag increases in the presence of polymers, contrary to what is found in single-phase flows. In other words, we show that the resulting drag is larger than what is observed in a Newtonian fluid with the same particle volume fraction . The particles have two main effects on the viscoelastic flow behavior: (i) Their presence limits the polymer stretching, hence the drag reducing effects of polymers are quenched, and (ii) the larger strain rates near the particle surface increase the particles and polymeric shear stresses. In addition, (iii) the shear-thinning behavior of the fluid enhances the particle migration towards the walls, which strengthens the near-wall particle layer and further increases the overall drag.
In the simulations, we consider an incompressible fluid governed by the Navier-Stokes equations, with the non-Newtonian feature modeled by the constitutive finite extensible nonlinear elasticity-Peterlin (FENE-P) closure [14], able to capture both the elasticity and shear-thinning behaviors of polymers. In particular, the fluid is governed by the incompressible Navier-Stokes equations with the additional viscoelastic stress τ E i j which accounts for the non-Newtonian nature of the flow, defined as τ E i j = [C i j /(1 − C kk /L 2 ) − δ i j ]/Wi. In the expression above, C i j is the conformation tensor, L the dumbbell extensibility, δ i j the Kronecker delta, and Wi the Weissenberg number. The rigid particle motion is described by the Newton-Euler equations, and their presence is modeled with an immersed boundary method (IBM) [15]. We use an explicit fractional-step method to solve the system of equations, where all the terms are advanced with the third-order Runge-Kutta scheme. We include global artificial diffusion in the transport equation of the conformation tensor in the first transient part of the simulation, and then gradually reduce it while approaching the final statistical steady state regime, when it is completely removed. Also, we use a fifth-order weighted essentially nonoscillatory (WENO) scheme for the advection terms in the conformation tensor equation. Apart from that, the governing differential equations are solved on a staggered grid using a second-order central finite-difference scheme. We use a Cartesian uniform mesh in a rectangular box of size 6h × 2h × 4h, discretized with 960 × 320 × 480 grid points, corresponding to 32 points per particle diameter. The flow is driven at a constant flow rate. The numerical code has been extensively used and validated in the past for particulate flows [16,17], turbulent Non-Newtonian flows [18,19], and multiphase problems in non-Newtonian fluids [20]. The rigid particles are assumed to be neutrally buoyant and spherical with radius R equal to one tenth of the channel half height h, i.e., R = h/10. Fixing the particle density and size, the flow is governed by five additional nondimensional parameters: the Reynolds number Re, the Weissenberg number Wi, the particle volume fraction , the viscosity ratio β, and the dumbbell maximum extensibility L. We define the Reynolds number Re = U b h/ν based on the bulk velocity U b , channel half height h, total kinematic viscosity ν, and the Weissenberg number Wi as the ratio between the elastic and viscous forces, i.e., Wi = λU b /h, with λ the polymer relaxation time. We performed simulations at a fixed bulk Reynolds number equal to 2800 with different Weissenberg numbers in the low drag reduction regime (LDR) [9,10] (Wi ∈ [0, 4]) and particle concentrations ( ∈ [0, 0.2]) such that we cover different cases: a Newtonian single-phase flow (Wi = 0 and = 0), Newtonian particulate flows (Wi = 0 and > 0), non-Newtonian single-phase flows (Wi > 0 and = 0), and non-Newtonian particulate flows (Wi > 0 and > 0). If not otherwise stated, L = 60 and β = 0.9.
Using the previous set of parameters, the single-phase Newtonian case has an average friction Reynolds number Re τ = u τ h/ν (with u τ the friction velocity) equal to 180 [21], while the polymer addition provides a drag reduction DR of approximately 9% for Wi = 2 and 18% for Wi = 4. Here, the drag reduction is defined as DR = (τ w N − τ w )/τ w N , where τ w is the total mean shear stress at the wall and τ w N the value for the Newtonian case at the same volume fraction. When we introduce rigid particles in the flow, the friction Reynolds number of the suspension grows monotonically for all the Weissenberg numbers Wi considered, as shown in Fig. 1. However, the friction Reynolds number in the non-Newtonian cases (Wi > 0) grows with the volume fraction faster than in the Newtonian fluid (Wi = 0), which results in drag reduction (DR > 0) for small volume fractions and a drag increase (DR < 0) for high particle concentrations. This result is supported by the recent experimental measurements in Ref. [22], reported in the inset of Fig. 1, obtained in similar conditions but in a square duct and at a slightly higher Reynolds number.
Frictional Reynolds number Re τ as a function of the particle volume fraction for different Weissenberg numbers Wi. The inset shows the corresponding drag reduction DR computed with respect to the Newtonian multiphase case with the same volume fraction. The black crosses in the inset are experimental data taken from Ref. [22]. In the figure, all viscoelastic cases are obtained with L = 60 and β = 0.9, except for selected cases with L = 190 and β = 0.99 ( ) and L = 600 and β = 0.999 ( ).
The drag reducing effect of polymers is due to their elastic behavior, in terms of stretching and elongation: Polymers are stretched by the fluctuating shear and dissipate the accumulated energy when relaxing their length back to equilibrium [23]. In the presence of semidense particle suspensions, this mechanism is inhibited by the motion of particles that hinders the capability of the polymeric chains to fully stretch. To prove this, we report in Fig. 2 the square root of the mean of the trace of the conformation tensor C i j , i.e., C ii , where the overline is used to indicate the mean value, normalized by the dumbbell maximum extensibility L. In the single-phase flow the polymers tend to strongly stretch in the near-wall regions, where we see the peak of the distribution of the  trace of the conformation tensor, while they remain essentially coiled in the bulk of the channel, where the distribution smoothly and monotonically decreases to the minimum at the centerline. The behavior in the presence of particles is completely different: Although in the very near-wall region a moderate elongation is still present (around 20% of the single-phase flow), across most of the channel height we observe an almost flat distribution with low mean values. This behavior is due to the presence of the particles, as also visually confirmed in Fig. 3. Indeed, in the vast majority of the channel the trace of the conformation tensor is small and roughly uniform also instantaneously, except in the very near-wall regions.
While the diagonal components of the conformation tensor reduce, the polymer shear stress increases. To confirm this, we examine the mean shear stress balance across the channel. For a plane channel flow, the mean shear stress τ 12 can be decomposed into the sum of the viscous shear stress τ V 12 = μdu/dy, the Reynolds shear stress τ R 12 = −ρu v , the viscoelastic shear stress τ E 12 , and the particle contribution τ P 12 , with their sum being a linear function of the distance from the wall, equal to the total shear stress at the wall and null at the center of the channel. In the previous relations, the prime symbols are used to distinguish the fluctuation with respect to the mean (overline). These shear stress contributions are shown in Fig. 4 for four representative cases, all normalized by the Newtonian single-phase wall shear stress τ w N . In the Newtonian case, we observe the classical  behavior, with τ V 12 being maximum at the wall and rapidly vanishing in the bulk of the channel and the Reynolds shear stress τ R 12 null at the wall and maximum in the near-wall region. The behavior is modified by the presence of the rigid particles; in particular, the total shear stress grows because of the particle stress contribution τ P 12 which is almost uniform across the channel, except for being null at the wall and exhibiting a maximum located at around y ≈ R. Moreover, the viscous stress τ V 12 increases at the wall while the Reynolds shear stress τ R 12 decreases in the bulk of the channel. When polymers are added in a particulate flow, the total shear stress further increases. In addition to the relatively small shear stress contribution from the polymers τ E 12 , the overall drag increase is related to the growth of the Reynolds shear stress τ R 12 and to the increase of the particle τ P 12 and viscous stress τ V 12 , especially in the near-wall regions. This increase of near-wall activity is due to the shear-thinning behavior of the fluid, as will be explained later.
The loss of the drag reducing effect is explained by the increase of the Reynolds stresses upon the addition of particles, which induce local instantaneous high shear rates in their surroundings [24,25]. We show this in Fig. 5 by the probability P of the local shear rateγ l for the same four cases discussed in Fig. 4. The presence of particles strongly increases the spectrum of shear rates by enhancing the probability of both weak and strong events, leading towards a more flat distribution. The presence of polymers in the particulate case further modifies the distribution by reducing the probability of low shear rates and increasing that of high shear rates. The increase of the local shear rate is not accompanied by substantial polymer stretching, as previously observed in Fig. 3, because of the different nature of the fluctuations in the suspension. It is known that in single-phase polymeric flows, turbulent fluctuations are overall damped and the flow is highly anisotropic, with strong streamwise velocity fluctuations u and reduced wall-normal ones v , as shown in the top panels of Fig. 6. This is associated with the strengthening of the elongated streamwise correlated streaks. On the other hand, the presence of particles induces an opposite effect: Streaks are disrupted by the particles, resulting in a decrease of u and an increase of v , overall leading to a more isotropic flow. The same happens in the multiphase non-Newtonian flow, as shown in the bottom panels of Fig. 6. Fluctuations and local shear rates are enhanced in the particulate non-Newtonian flow mostly in the homogeneous directions, where the mean shear is null, thus not leading to substantial polymer stretching. Not only the fluid but also the particle dynamics are modified by the presence of polymers; in particular, we measured a 7% increase of the particle kinetic energy in the presence of polymers, which is consistent with the larger particle shear stress previously observed, confirming a picture where the coupling of polymers and particles is detrimental in terms of drag reduction due to an increased particle interaction.
As seen, the presence of the particles in a viscoelastic fluid reduces the polymer stretching and the consequent drag reduction and induces an increase of the total shear stress also by the growth of the polymeric shear stress, due to regions of high strain close to the particle surface. One additional effect responsible for further drag increase arises when the fluid is shear thinning: particle migration towards the wall. We show the local mean particle concentration l across the channel in Fig. 7  and observe that particles tend to concentrate in the middle of the channel [17] and in a nearwall layer, as demonstrated by the peaks of the local volume fraction at y = h and y ≈ R [16]. In particular, the latter expresses the tendency of finite-size rigid particles to accumulate close to the walls, which has been recognized to be mainly responsible for the increase in drag in turbulent suspension flows. In the viscoelastic cases, the concentration in the bulk of the channel reduces while the near-wall peak increases, inducing a further increase in the overall drag; this behavior is caused by the shear-thinning property of the fluid, which is well known to induce particle migration towards the wall [26]. Indeed, in viscoelastic Poiseuille flows two different attractors exist for a suspended particle, one at the wall and one at the centerline. The former depends on the level of shear thinning with pronounced shear thinning increasing the wall attraction. To quantify this effect on the drag modifications, we have performed additional simulations at = 0 and 0.2 with different levels of shear thinning, displayed by the rheological curves in the inset of Fig. 7. In particular, we have varied β and L keeping the product L 2 (1 − β ) constant in order to have a similar amount of drag reduction as in the single-phase case [27]. As reported in Fig. 1, we observe that the drag decreases as the shear-thinning effect is reduced in the particulate case at = 0.2, but remaining larger than the one measured in a Newtonian solvent in the range of parameters studied here. The reduced level of drag increase is due to the loss of the particle migration: Indeed, in Fig. 7 we can observe that the particle concentration profile in the viscoelastic case with β = 0.999 is approximately equal to what is found in the Newtonian carrier fluid.
To conclude, we have studied viscoelastic turbulent flows of particle suspensions by means of direct numerical simulations. We have shown that, in the presence of a moderate concentration of particles in the low drag reduction regime, the beneficial effect of polymers in terms of drag reduction are lost, and an increase in drag is observed for sufficiently dense suspensions. The inhibition of the drag reducing mechanism in the presence of particles is due to the quenching of the stretching of the polymers which essentially nullifies the elastic behaviors of the fluid. On the other hand, the further increase in drag is due to an increase of the particle and polymeric shear stresses, especially close to the particle surface, and to the strengthening of the particle near-wall layer induced by their migration towards the wall due to the shear-thinning fluid behavior. The combination of these three effects explains the increase of drag in semidense and dense particulate suspensions in viscoelastic fluids. L.B. was supported by the ERC-2013-CoG-616186 TRITOS grant and by the Swedish Research Council, VR 2014-5001, and acknowledge the computer time provided by SNIC (Swedish National Infrastructure for Computing).

APPENDIX: CODE VALIDATION
The numerical code used in the present Rapid Communication has been extensively used and validated for single-phase and multiphase flows of shear-thinning, viscoelastic, and elastoviscoplastic fluids with more details on the algorithm and validation campaign discussed in previous publications [17][18][19][20]28,29].
For the sake of completeness, we provide here a further validation of our code by comparing the angular velocity of a sphere in an homogeneous viscoelastic simple shear flow, described by the Oldroyd-B model, against experimental and numerical data found in the literature [30,31]. We chose a rectangular domain of size 8R × 4R × 4R in the streamwise x, wall-normal y, and spanwise z directions, with R the particle radius. The numerical domain is discretized with 32 grid points per particle diameter as in the rest of the work. The particle is initially positioned at the center of the box and a constant shear rate (γ 0 ) is applied, generated by the two parallel walls moving in the x direction. Periodic boundary conditions are considered in the y and z directions. The particle Reynolds number Re is set equal to 0.025 and the Weissenberg number Wi is varied in the range from 0 up to 2. Figure 8 reports our results and those from the literature; as can be seen, there is an excellent agreement between our numerical results and the results previously obtained, thus confirming the validity of the present numerical implementation.