Simulations of a sphere sedimenting in a viscoelastic fluid with cross shear flow

Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60

Contents lists available at SciVerse ScienceDirect

Journal of Non-Newtonian Fluid Mechanics journal homepage: http://www.elsevier.com/locate/jnnfm

Simulations of a sphere sedimenting in a viscoelastic fluid with cross shear flow S. Padhy a,⇑, E.S.G. Shaqfeh a,b, G. Iaccarino a, J.F. Morris c, N. Tonmukayakul d a

Department of Mechanical Engineering, Stanford University, Stanford, CA 94305, USA Department of Chemical Engineering, Stanford University, Stanford, CA 94305, USA c Levich Institute and Department of Chemical Engineering, The City College of New York, New York, NY, USA d Halliburton, Duncan, OK 73536, USA b

a r t i c l e

i n f o

Article history: Received 24 September 2012 Received in revised form 7 February 2013 Accepted 8 February 2013 Available online 26 March 2013 Keywords: Viscoelasticity Settling rate Shear flow Particle sedimentation Suspension

a b s t r a c t The settling rate of heavy spheres in a shear flow of viscoelastic fluid is studied by numerical simulation. Experimental data [Tonmukayakul et al., US Patent Application US20110219856 (2011); van den Brule and Gheissary, J. Non-Newton. Fluid Mech. 49 (1993) 123–132] have shown that both shear thinning and the elasticity of the suspending polymeric solutions affect the settling rate of the solids. In the present work, simulations of viscoelastic flow past a single, torque-free sphere with a cross shear flow are used to study the effect of the elasticity of the carrying fluid on the sphere’s settling rate. The FENE-P constitutive model is used to represent a viscoelastic Boger fluid, with parameters obtained by fitting rheological data. A twofold increase in drag, i.e. a decrease in settling rate, is obtained with increase in the cross shear Weissenberg number, Wi 6 15, even though the shear viscosity of the solution decreases over this same range. At small Weissenberg number, Wi < 2, the simulations remain in quantitative agreement with the experiments. At higher Weissenberg number, the numerical results remain in qualitative agreement with settling experiments although the magnitude of the simulated decrease in settling rate is smaller than that observed in experiments. The detailed physical mechanism for the increase in the drag experienced by the sphere in the simulations is presented and we show that s11 component of the viscous stress (with 1, the sedimentation direction) is the primary cause of the increase in sphere drag. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Suspensions of rigid particles in viscoelastic fluid play an important role in many applications [1]. The example motivating this work is the hydraulic fracturing operation used in oil and gas wells where suspensions of solid particles in polymeric solutions are pumped into hydraulically-induced fractures. The particles which are usually dense mineral materials like sand or sintered bauxite must prop these channels open to enhance the rate of petroleum recovery [2]. In a vertically-oriented fracture, sedimentation of the particles thus limits the axial distance traveled by the so-called ‘proppant’ particles as they drop to the bottom of the fracture. Given the need to limit settling in order to place particles far into the fracture, one might expect that a liquid of high viscosity would be employed in the fracturing operation, but this is not an option. Because fractures are deep below the Earth’s surface, use of a ⇑ Corresponding author. Address: 488 Escondido Mall, Building 501, Room 501S, Stanford, CA 94305-3035, USA. Tel.: +1 6503538805; fax: +1 6507253525. E-mail addresses: [email protected] (S. Padhy), [email protected] (E.S.G. Shaqfeh), [email protected] (G. Iaccarino), [email protected] (J.F. Morris), [email protected] (N. Tonmukayakul). 0377-0257/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jnnfm.2013.02.003

high-viscosity liquid to limit settling leads to excessive pumping pressures to drive flow in the well-bore and in the fracture itself. In fact, a commonly used proppant-transport liquid is an aqueous guar polymer solution, transiently cross-linked with borate ion [3]. This relatively low-viscosity liquid exhibits significant shear thinning (viscosity reduction with increasing shear rate), yet has proven effective in transporting the dense solid particles into the fracture. Here, we examine how the non-Newtonian, and in particular the elastic, properties of the liquid provide this particle-transport capability. It is well-known that particle settling in a viscoelastic fluid can be quite different from that which is observed in Newtonian fluids. An issue of importance in this context is that in Newtonian fluids, the presence of an imposed shear flow in the direction perpendicular to gravity—which we term a cross shear flow—has no effect on the settling of a spherical particle in Stokes flow. By contrast, in a non-Newtonian liquid, the complex rheological properties induce a nonlinear coupling between the sedimentation and shear flow. It is striking that this problem has gone unexplained for nearly two decades. Van den Brule and Gheissary [4] performed experiments on single sphere settling in the vorticity direction of an imposed shear flow in aqueous mixtures containing polyacrylamide

49

S. Padhy et al. / Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60

to introduce elasticity, and found reductions of settling velocity as large as 50% simply due to modest shearing of the liquid. We present our own experimental results [5] to show that the settling rate of more concentrated suspensions in viscoelastic solutions can likewise be reduced by a cross shear flow (the reduction level surpasses that reported by van den Brule and Gheissary [4]), thus establishing that the phenomenon of reduced settling is not solely one of isolated particles. Viscoelastic reduction of settling rate in flowing slurries is thus of direct relevance to applications involving realistic loadings of solids. The influence of individual particle effects on rheology of suspensions in a polymeric solution has gained significant attention [6]. Both experimental [7] and numerical studies [8,9] have examined the effect of fluid viscoelasticity on the rotation of a nonBrownian sphere in simple shear flow (x), showing its reduction from the value in a Newtonian fluid in Stokes flow, i.e. x ¼ c_ =2 where c_ is the bulk shear rate. An analytical result for the angular velocity derived by Housiadas and Tanner [10] is in agreement with experimental and numerical studies. Goyal and Derksen [11] also verified this result using direct numerical simulation based on lattice-Boltzmann scheme. The behavior of multiple spherical particles in a viscoelastic medium under shear flow, without gravitational settling, has been studied with a focus on the microstructure [12–14]. The particles were found to form string-like structures aligned along the shear flow direction. An asymptotic theory for viscoelastic fluids [15] predicts the chaining of weakly interacting particles, both for spheres sedimenting without shear flow as well as spheres in shear flow in the absence of gravity. There have also been studies of sedimentation in viscoelastic fluids. Bobroff and Phillips [16] used nuclear magnetic resonance imaging to study sedimentation rates of systems of spherical particles in various non-Newtonian fluids. They found that sedimentation rate becomes slower with time, while elongated columns of particles form in the direction of gravity. Similar non-homogeneous microstructures have been observed in other experiments [17,18]. Goyal and Derksen [11] performed numerical simulations of sedimentation of single particle as well as multiple particles in viscoelastic fluid. They found the interaction between spheres depends strongly on the elastic properties of the liquid. There has, however, been little work devoted to the issue of interest here, namely sedimentation under imposed shear for viscoelastic fluids. Recent experiments [5] have shown a reduction of settling rate by shear flow for assemblies of spherical particles in various viscoelastic liquids. Similar to the findings of van den Brule and Gheissary [4] for a single sphere, this effect is observed for Boger fluids, where viscosity is independent of shear rate, and for shear-thinning liquids. Housiadas and Tanner [19] examined this problem for a weakly viscoelastic fluid using perturbation theory. Their work provides some insight on the effect of elasticity in a weakly viscoelastic fluid on sedimentation. There has not been any definitive quantitative explanation of the mechanisms by which elasticity in the fluid affects sedimentation. We consider the work by Housiadas and Tanner [19] in some detail in Section 5.2 and following. In the present study, we perform numerical simulation of the fluid motion and compute the stresses caused by settling of a spherical particle in a cross shear flow in order to develop a mechanistic understanding of how viscoelasticity reduces the sedimentation rate of dense particles. We model the rheology of fluids from the studies of van den Brule and Gheissary [4] and from Tonmukayakul et al. [5] to obtain the parameters for the simulations. The sphere drag, computed from the simulations, is directly related to the settling rate of the sphere. The numerical results indicate that reductions of the sedimentation speed of up to an order of magnitude may be achieved for sufficiently large (but potentially experimentally accessible) Weissenberg number, Wi, a measure of the elasticity of the fluid.

2. Problem formulation 2.1. Problem definition We consider a freely sedimenting sphere moving with constant speed U1 under the action of gravity in a cross shear flow. Simulations for this problem are performed in an inertial frame of reference attached to the sphere. Thus, in this frame of reference, a ~i is a unit fluid with uniform velocity of Ui (where U i ¼ U 1 g~i ; g vector in the direction of gravity and from Fig. 1, g~i  di1 ) is flowing past a stationary sphere. The cross shear flow is applied by enforcing a no slip boundary condition at the upper and lower plates and moving them with a speed S in opposite directions. A rotational velocity, xi, is applied to the sphere subject to the constraint that the torque on the sphere is zero. By symmetry, only x1 is nonzero. 2.2. Governing equations The dimensionless governing equations for the flow of an incompressible fluid containing polymer additives are obtained by conservation of mass and momentum, viz.

@ui ¼0 @xi

ð1Þ

@ui @ui @p b @ ui 1  b a @s ¼ þ þ þ uj @xi Re @xj @xj Re Wi @xj @t @xj 2

P ij

ð2Þ

where sPij is the stress due to the polymers. The velocity scale, U1, is taken as the sedimentation speed of the sphere. The shear Weissenberg number, Wi ¼ kc_ , is defined as the product of the characteristic relaxation timescale of the polymer k and the characteristic shear rate c_ applied across the cell in the x2  x3 plane. The Reynolds number, Re = qU1D/(ls + lp), is defined as the ratio of the inertial forces to viscous forces where D is the sphere diameter, ls is the solvent contribution to the viscosity and lp is the polymer contribution to the viscosity with b = ls/(ls + lp). Another dimensionless parameter of interest is the ratio of the characteristic shear rate, c_ and U1/D, thus a ¼ c_ D=U 1 , which can be combined with Wi to define the flow Weissenberg number, h = kU1/D = Wi/a. Note h is defined as the ratio of the characteristic relaxation timescale and the characteristic flow timescale D/U1. The FENE-P model is used as the constitutive equation. The polymer stress term in Eq. (2) is given by [20],

sPij ¼

cij  dij 1  cLkk2

ð3Þ

where L is the dimensionless maximum polymer extensibility and cij is the dimensionless polymer conformation tensor with both scaled by the equilibrium Hookean spring length. The conformation tensor is governed by the transport equation given by,

@cij @cij @uj @ui a P 1 @ 2 cij  cik  ckj ¼ s þ þ uk @t @xk @xk @xk Wi ij ReSc @xk @xk

ð4Þ

where the Schmidt number, Sc = (ls + lp)/(qC), is defined as the ratio of momentum diffusivity to spatial diffusivity denoted by C.

Fig. 1. Schematic of the problem used for simulation.

50

S. Padhy et al. / Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60

Note that a constant artificial dissipation for the cij equations may be introduced for numerical stability by adding an extra diffusion term to the right-hand side of Eq. (4). This extra diffusion term smears the sharp gradients of polymer stress that can form due to the hyperbolic nature of Eq. (4). Many studies have used similar artificial global stress diffusion with Sc 6 1 [21] and Sc  10 [22] for numerical stability. We demonstrate the effect of this diffusion term on our calculations in Sections 5.1 and 3.1. We show that the changes in drag coefficients are minimal for Sc P 106 for Re = 0.01. Our simulations results are presented for Sc = 1 barring a few simulations which have Re = 0.01, Sc = 106. The single mode Giesekus model [23] is used as the constitutive equation for the comparison of our results for the rotation rate of a sphere in a viscoelastic shear flow with Snijkers et al. [9] (Section 3.3). In this case, the polymer stress term in Eq. (2) is given by,

The conformation tensor is governed by the transport equation given by,

Since, as described above, a small amount of artificial diffusion may be added for numerical stability, in which case the boundary conditions for cij need to be explicitly specified at the no slip boundaries for the simulations with the artificial diffusion term. A no flux condition is specified in the normal direction for cij for the top and bottom plate as well as at the sphere (@cij/@n = 0). We have demonstrated in previous publications [22,24] this boundary condition makes small difference for Sc  1. A convective outlet condition is applied for cij at the outlet. At the inlet, the analytical solution of cij for plane Couette flow of a FENE-P fluid is imposed [27]. In the x3 direction a periodic boundary condition is applied both for velocity and cij. All the simulations were carried out at Re = 0.01. The time step size Dt was chosen to be 4.261  105k. The results we discuss below represent the steady state results. The steady state condition was determined to be reached when the drag value on the sphere as well as the L1 norms of s11, sP11 ; sP12 and sP13 in the domain changed by less than 1% in 23,000 1 k time steps (k). Here s1k ¼ @u þ @u is the viscous stress on the @xk @x1 sphere.

@cij @cij @uj @ui a P af P P 1  cik  ckj ¼ s  s s þ þ uk @t @xk @xk @xk Wi ij Wi ik kj ReSc

2.4. Mesh

sPij ¼ cij  dij

ð5Þ

@ 2 cij  @xk @xk

ð6Þ

where f is the dimensionless mobility parameter [23]. The Eqs. (1) and (2) along with either (4) or (6) provide a set of 10 coupled partial differential equations for obtaining the instantaneous velocity, pressure, polymer conformation and stress field. The numerical method we use for solving the equations is described in detail in Richter et al. [22]. The viscoelastic code is built upon an existing incompressible Newtonian flow solver developed at Stanford’s Centre for Turbulence Research named CDP. It is a fully three dimensional code based on an unstructured finite volume formulation and has been used efficiently on a variety of parallel architectures. The polymer stress term (either FENE-P or Giesekus) is added in the Newtonian momentum solver as a volumetric source term and the conformation evolution equations are solved as six coupled scalar transport equations. This code is capable of solving problems with a wide range of Reynolds and Weissenberg numbers [22,24–26]. 2.3. Boundary conditions The standard domain (with directions as defined in Fig. 1) has dimensions of 37.5D, 11.66D and 5.83D in the x3, x2 and x1 directions, respectively. The sphere is placed at the origin of the coordinate system. The boundary where fluid is entering the domain in the x1 direction has a velocity of di1  (ax2)di3. A no slip boundary condition is applied at the top and bottom plate in the x2 direction with velocities of the plates given by di1  5.83adi3 and di1 + 5.83adi3 respectively. A convective outlet condition is applied where the fluid is leaving the domain in the x1 direction. The velocity there is extrapolated from the interior cells assuming fully developed flow and maintaining the overall mass balance. A fixed velocity is imposed at the surface of the sphere corresponding to the solid body rotation of the sphere with rotational velocity x1 so that the net torque on the sphere is zero. The value of x1 is found by iteration using the secant method until the magnitude of the non-dimensional torque value ðT ¼ 2T=ðpðls þ lp ÞD3 c_ Þ, where T is the dimensional torque value) becomes less than 0.01. The rotational velocity is updated after every 23,000 time steps (k) according to the Eq. (7) where k is the time step number.

x1ðkþ1Þ ¼ x1ðkÞ 

ðk1Þ xðkÞ 1  x1

T ðkÞ  T ðk1Þ

T ðkÞ

ð7Þ

Simulations were carried out on two different unstructured meshes to check grid convergence and the effect of the domain size. One of the mesh has dimensions of 37.5D, 11.66D and 5.83D in the x3, x2 and x1 directions, respectively and contains 6.268 million tetrahedral elements. The other mesh has dimensions of 37.5D, 12.5D and 11.66D in the x3, x2 and x1 directions, respectively and contains 5.34 million hexahedral elements. Both meshes have smaller elements near the sphere to resolve the polymer stress gradients and larger elements farther from the sphere. The tetrahedral and hexahedral meshes have a smallest element size of 0.0125D and 0.0210D respectively. The results from the simulations carried out in both the meshes are in agreement (Section 5.1). The mesh with the tetrahedral elements is used for all the simulation results for Wi < 10 presented in the paper unless otherwise specified. For higher Weissenberg numbers, a mesh with larger domain size is needed. The mesh that is used for Wi P 10 has dimensions of 37.5D, 11.66D and 10D in the x3, x2 and x1 directions, respectively and contains 9.294 million tetrahedral elements. It has the smallest element size of 0.01D near the surface of the sphere.

3. Verification studies The verification of the code is carried out by comparison with previous results in the literature for similar work. We begin with comparison of our simulations for the steady, fully-developed Poiseuille flow in a two-dimensional channel including a known analytical result [28]. We proceed to the comparison of the drag coefficient of sedimentation of a sphere along the centerline of a tube [29] and the rotation rate of a sphere in a viscoelastic shear flow [9]. It is important to note that in these problems both the shear Weissenberg number (Wi) as well as the uniform flow Weissenberg number (h) play an important role. 3.1. Poiseuille Flow in 2D channel An analytical solution for the steady, fully-developed Poiseuille flow of viscoelastic fluid with FENE-P model has been previously derived by Cruz et al. [28]. This result was used by Richter et al. [22] for code verification. Our simulations for flow of viscoelastic fluid in a channel is carried out in a mesh with quadrilateral elements of size 0.01h where h is the separation between the two

51

S. Padhy et al. / Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60

infinite parallel plates. The parameters for the FENE-P model are chosen to be L = 100 and b = 0.9 as a test case. The average velocity magnitude ðUÞ and channel half-width (h/ 2) are used to make the variables non-dimensional. Thus, the flow Weissenberg number is given by h ¼ 2kU=h and the Reynolds number is given by Re ¼ qUh=ð2ls þ 2lp Þ. For the simulations, we first use parameters h = 10, Re = 0.01 and Sc = 1 (i.e. spatial diffusivity, C = 0). The comparison between the computed and analytical results for velocity and conformation tensor components is shown in Fig. 2a–c. The computed results match the analytical profiles very well. It is important to note that addition of artificial diffusion term in the Eq. (4) is not necessary for convergence of results in this case. Thus, this case is also used to test the effect on accuracy of the solution with the addition of artificial diffusion and we now present results over the same parameters but varying Sc. For the comparison, the relative error() between computed and exact solution is computed using:

  kYsim k2  kYexact k2    kYexact k2

 ¼ 

ð8Þ

3.2. Sedimentation of a sphere in a viscoelastic fluid Yang and Khomami [29] studied the problem of determining the drag coefficients for a sphere sedimenting in a tube using different numerical models. They compared their results with the experimental results obtained by Arigo et al. [30]. We chose one of their computational set of results for comparison to verify our code. The comparison is carried out for a sphere-to-tube radius ratio of a/R = 0.243 with the FENE-P parameters of g = lp + ls = 13.76 Pa s, k = 0.7925 s, L = 30 and b = 0.59 given in Arigo et al. [29]. The geometry used for this comparison is a sphere inside a cylindrical tube where the length of the tube in the direction of sedimentation is 30a. The sphere is placed at a distance of 4a from the inflow boundary. The mesh contains 1 million hexahedral elements with smaller elements near the surface of the sphere. The non-dimensional parameters are calculated using the sedimentation speed of the sphere, U1, as the velocity scale and radius of the sphere, a, as the length scale. Thus, the characteristic strain rate is given by U1/a and the flow Weissenberg number is given by h = kU1/a. The Reynolds number is given by Re = qU1a/(ls + lp). The shear Weissenberg number, Wi =0, as there is no cross shear flow in this problem. The simulations are carried out with Re = 0.05 and Sc = 1. The viscoelastic drag coefficient, KVE, is calculated using:

1

1

0.5

0.5

0

0

y

y

where kk2 denotes the L2 norm of the variable and Y denotes the conformation tensor components. The relative errors in c11 and c12 components of the conformation tensor is calculated over a large range of Schmidt and Reynolds number. Other parameters are kept fixed i.e. h = 10, Wi = 0, L = 100 and b = 0.9. The relative error for c11 and c12 as a function of ReSc is shown in Fig. 2d. The figure shows results for two different Reynolds number and various Schmidt numbers. The relative errors for different Re are nearly equal once they are plotted as functions of ReSc instead of just Sc. Thus, for evaluating the effect of artificial diffusion ReSc is the important parameter and not just Sc. This is also confirmed by the fact that ReSc appears in the denominator of the artificial diffusion term in

Eq. (4). The relative errors decrease as ReSc increases as the errors introduced due to the artificial diffusion term is reduced. At ReSc = 1, the error due to artificial diffusion term is zero and the relative error is now determined by the numerical diffusion due to the discretization. The change in relative error is very small as ReSc changes from 104 to 1. Thus, the effect of artificial diffusion is negligible and the error is dominated by numerical diffusion due to discretization for ReSc P 104.

−0.5 −1

−0.5

0

2

4

6

−1

8

(a)

(b)

4000

6000

1Ε6

∞

−1

Relative Error (c11,c12)

10

0.5

y

2000

c11

1

0 −0.5 −1 −40

0

u (y)

−2

10

−3

10

−4

−20

0

c12

(c)

20

40

10

1Ε2

1Ε4

ReSc

(d)

Fig. 2. Comparisons of analytic (—) and simulated () (a) streamwise velocity u and conformation tensor components, (b) c11, and (c) c12 for Re = 0.01, Sc = 1, h = 10, b = 0.9 and L = 100. (Note that the simulated result for only 1 in every 5 mesh elements is plotted for clarity in the plot.) (d) Relative Error in c11 and c12 for different Re. , c11Re = 0.01; , c11 Re = 1; , c12 Re = 0.01; , c12 Re = 1.

52

K VE ¼

S. Padhy et al. / Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60

F1 ðlp þ ls ÞU 1 a

ð9Þ

where F1 is the force experienced by the sphere in the x1 direction. The Newtonian drag coefficient, KN, is calculated using Faxen’s law. It is given by 6pc where the correction factor c is 1.93 [31]. The result for the normalized drag coefficient (the ratio of viscoelastic drag to Newtonian drag) as a function of Weissenberg number is shown in Fig. 3a. The agreement between the simulation results obtained by a finite element code by Yang and Khomami [29] and our finite volume code is found to be quite reasonable. Our simulation results also shows the initial dip in the normalized drag coefficient before increasing with Wi which is consistent with the results obtained by Yang and Khomami [29]. We also carried out the simulations on a different grid and at Re = 0.01.The results stayed nearly the same and at this point the small difference in the simulation results can only be attributed to the difference in mesh resolutions. 3.3. Rotation of a sphere in viscoelastic shear flow Snijkers et al. [9] studied the problem of determining the rotation rate of a sphere in a viscoelastic shear flow using numerical simulations and experiments. They obtained the rotation rates in a Boger fluid experimentally and used the single mode Giesekus model for obtaining the simulation results. The parameters for the model were determined by fitting to the rheological data. We used the same parameter values for our simulations as determined by Snijkers et al. [9] (g = lp + ls = 45.8 Pa s, k = 6.18 s, f = 0.00032 and b = 0.611). The velocity scale in this problem is defined as the magnitude of the velocity with which the top and bottom plate translates in the x3 direction. This is denoted by S as shown in the Fig. 1. The length scale is defined as the diameter of the sphere (D). The characteristic shear rate is given by c_ ¼ 2S=H where H is the separation between the plates in the x2 direction. Thus, the shear Weissenberg number is given by Wi ¼ kc_ and the Reynolds number is given by Re = qSD/ (ls + lp). The flow Weissenberg number, h = 0, as there is no uniform flow in this problem. The simulations are carried out with Re = 1 and Sc = 106. The comparison between our results and those obtained by Snijkers et al. [9] is shown in Fig. 3b. The agreement is very good. The maximum error between the simulation results by Snijkers et al. [9] and our simulation results is 2.6%. The agreement between the simulation and experimental results is especially good at low Weissenberg number. The rotation rates predicted by simulations are slightly higher than experimental values for Wi P 4. This might be due to the fact that single mode Giesekus model is not sufficient to model the rheology for Wi P 4.

4. Fitting FENE-P model parameters to experimental data The parameters for the FENE-P model were obtained by fitting the experimental rheological properties of 1% PAA solution (in a mixture of water and corn syrup) measured by Tonmukayakul et al. [5] using an L2 error measure. The relaxation time is obtained by fitting the experimental data for first normal stress difference obtained by the stress relaxation test. The parameters so obtained after the fitting procedure were g = lp + ls = 1.249 Pa s, k = 0.469 s, L = 31.6 and b = 0.785. The comparison of the fit and experimental data is shown in Fig. 4. Note that the results correspond to modest shear thinning of g over the experimental range of shear rate. A similar procedure is followed for 100 ppm PAA solution in glucose syrup reported by van den Brule and Gheissary [4]. The parameters obtained were g = 5.94 Pa s, k = 0.284 s, L = 5.47 and b = 0.69. The comparison of the fit and experimental data is shown in Fig. 5. 5. Results and discussion Tonmukayakul et al. [5] measured the settling rate for a system of spherical particles at solid volume fraction of c = 0.2 in a viscoelastic medium under the action of cross shear flow. Van den Brule and Gheissary [4] measured the same for a single sphere. Both found a decrease in settling rate as shear rate increased. This trend has also been obtained by Housiadas and Tanner [19] analytically for weakly viscoelastic fluids. The settling rate is inversely proportional to the drag experienced by the particles. In the simulations, we compute the drag for a fixed velocity. As the settling rate is inversely proportional to the drag, the decrease in settling rate is equivalent to an increase in the drag coefficients in the simulations. The drag coefficients are calculated as the cross shear rate, and consequently the Weissenberg number, is increased. The drag coefficient is defined by:

Cd ¼

2F 1 ðlp þ ls ÞU 1 D

ð10Þ

where F1 is the force experienced by the sphere in the x1 direction. The simulation results are verified to converge to the correct solutions through various tests. The results for the convergence tests are presented in the next section. 5.1. Convergence results The convergence to steady state for h = 0.1 and Wi = 2.34 is shown in Fig. 6a. The temporal variation of L1 norms of s11 and sP13 in the domain is shown. The figure shows the temporal plot for Sc = 106 and Sc = 1 (no artificial diffusion term). The initial

1.3

1.5

1.2

ωλ

KVE/KN

1 1.1

0.5 1 0.9

0

1

2

3

4

5

0

0

2

4

θ

Wi

(a)

(b)

6

Fig. 3. Comparison to previous results (a) normalized drag coefficient for sedimentation of a sphere. a/R = 0.243, L = 30. Solid line, simulation results of Yang and Khomami [29]; }, our simulation results; , experimental results of Arigo et al. [30]. (b) Rotation rates of a sphere in viscoelastic shear flow. Solid line, simulation results of Snijkers et al. [9]; , our simulation results;h, experimental results for Boger fluids [9].

53

S. Padhy et al. / Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60 2

1.5

10

1

10

N1 (Pa)

η (Pa.s)

1.4 1.3 1.2

0

10

−1

10

1.1

−2

1 10

0

10

10

1

0

1

2

Wi

Time (sec)

(a)

(b)

3

Fig. 4. Fits to experimental data of 1% PAA solution (a) steady state viscosity and (b) transient first normal stress difference. Solid line, FENE-P model; j, experimental data of Tonmukayakul et al. [5].

4

10

6

3

N1 (Pa)

η (Pa.s)

5.5

5

10

2

10 4.5

1

0

10

1

10

0

1

10

10

10

Wi

Wi

(a)

(b)

Fig. 5. Fits to experimental data of 100 ppm PAA solution (a) steady state viscosity. Solid line, FENE-P model; j, experimental data of van den Brule and Gheissary [4]. (b) Steady state first normal stress difference. solid line, FENE-P model; j, experimental data of van den Brule and Gheissary [4].

Table 1 Convergence to steady state for h = 0.1, Wi = 4.69. Time

% Change in L1 norm of s11

% Change in L1 norm of

9k 10k 11k 12k

0.2 0.18 0.08 0.007

3.61 2.03 1.38 0.9

sP13

50

50

40

45

30

40

Cd

L Norm of Stresses

conditions for cij and the velocity are the analytical solution for plane Couette flow of a FENE-P fluid. The results reach steady state after approximately 10k. There is negligible difference between the temporal plots for both Schmidt numbers. The time convergence

for h = 0.1 and Wi = 4.69 is shown in Table 1). The simulation is terminated when the norms for all the stresses change by less than 1% which requires 12k for this case. Simulations were carried out using both the mesh with hexahedral elements (coarse mesh) and the mesh with tetrahedral elements (fine mesh). The results are shown in Fig. 6b. It can be seen that grid convergence has been achieved and the finer mesh is used for all further simulations. The effect of introduction of the artificial diffusion is also shown in Fig. 6b. The variation of drag coefficients with ReSc is shown. It is important to note that the effect of artificial diffusion in channel flow (Section 3.1) is negligible for ReSc P 104. This gives us an estimate for the magnitude of the Schmidt number to be chosen such that the effect of artificial

20

35

10

30

0

0

1

2

3

4

25 1Ε0

1Ε2

1Ε4

time

ReSc

(a)

(b)

∞

Fig. 6. (a) Convergence to steady state result for h = 0.1, Wi = 2.34. ——, s11 (Sc = 106); h, s11 (Sc = 1); - - -, sP13 ðSc ¼ 106 Þ; ; sP13 ðSc ¼ 1Þ. Grid convergence and convergence with ReSc, h = 0.1. (b) Wi = 4.69. , Mesh1 (hexahedral elements); , Mesh2(tetrahedral elements).

54

S. Padhy et al. / Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60

diffusion is minimal. As the simulations for this problem are carried out at Re = 0.01 thus a Schmidt number of 106 is sufficient to minimize the effect of artificial diffusion. This is further verified from the results shown in Fig. 6b. The changes in drag coefficients are 1.7% and 0.21% as Sc is increased from 105 to 106 and from 106 to 1 respectively. Thus, the changes in drag coefficients for Sc P 106 (at Re = 0.01) are negligible. 5.2. Drag coefficient results The drag coefficients depend on the shear Weissenberg number (Wi) as well as the uniform flow Weissenberg number (h = kU1/D). The simulations for the comparison with the experiments are carried out at the same values of h and Wi as the experiments. It is important to note that all the simulations were performed without the artificial diffusion term (Sc = 1). Some of the results for the simulations with Sc = 106 are also presented and compared to the results with Sc = 1 to show that the effect of artificial diffusion is negligible for Sc P 106. In Fig. 7, the comparison between our simulation results and the experimental results by van den Brule and Gheissary [4] is presented. In the figure, Cd/C0, is plotted as a function of Weissenberg number where C0 is the drag coefficient at Wi = 0 and h = 0.418. The simulations for this comparison were performed in a domain with 5D as the dimension in x2 direction to match the experimental conditions. The drag coefficients obtained from the simulations increase with Weissenberg number i.e. the settling rate decreases with Wi, as observed in the experiments. The simulation results show good agreement with the experimental results in the range of Weissenberg number (Wi < 2) used by van den Brule and Gheissary [4]. In Fig. 7, the comparison between our simulation results and the theoretical result for h 1 by Housiadas and Tanner [19] is also presented. The simulations for this comparison were performed in a larger domain with 11.66D as the dimension in x2 direction, as the theory is only valid for an infinite domain. We used the Oldroyd B model and matched the value of the parameters in the simulations to that used in [19] (b = 0.7 and k = 0.3 s). The simulation results are in good agreement to the theoretical results at small Weissenberg numbers, deviating for Wi P 0.6. This is expected as the theory is valid for small Weissenberg numbers only. The theoretical drag coefficients are also significantly smaller than the experimental results. We speculate at this point that this difference is primarily due to the effect of the walls on the sphere which is not taken into account in the theory. This is also suggested by the fact that the simulation results performed in the smaller

domain matches the experimental results. Further work is required to investigate how the walls affect the drag coefficient of the sphere. The experiments by Tonmukayakul et al. [5] were done at higher Weissenberg numbers. The results for higher Weissenberg numbers (Wi P 2) are shown in Fig. 8. In Fig. 8a and b, the ratio of drag coefficient, Cd/C0, is plotted as a function of Weissenberg number. C0 is the drag coefficient at Wi = 0, h = 0.1 for Fig. 8a and at Wi = 0, h = 0.328 for Fig. 8b (the experimental value of h is 0.328 for Wi = 0). For both our simulations and the experiments, there is an increase in drag coefficient (decrease in settling rate) as the Weissenberg number, Wi, is increased. However the increase in the drag predicted in the simulations is smaller than that observed in the experiments. This difference might be anticipated as the volume fraction used in the experiments (c = 0.2, Tonmukayakul et al. [5]) is much larger than that used in the simulations. In Fig. 8b, the simulation results without the artificial diffusion term is also presented. There are negligible changes in the results as Schmidt number is increased from 106 to 1. The mechanism for the increase in drag is further examined through simulations of a single sphere with h = 0.1. The drag coefficient contributions are broken down into the individual components associated with the stress contributions on the sphere. The total drag coefficient is given by:

 Z Z Z Z  @u1 @uk nk dS C d ¼ 2Re pn1 dS þ 4b þ @xk @x1 S S |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl ffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} form drag

Z Z 1b þ2 a sP1k nk dS Wi S |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

where nj is the unit outward normal vector on the surface of the sphere, p is made dimensionless with qU 21 and dS is the non-dimensional differential area element (non-dimensionalized by D2). The individual contributions to the drag with respect to the contributions at Wi = 0 are presented in Fig. 8c. A key finding is that the drag contribution due to the viscous stresses increases with Wi. The drag contribution due to the polymer stresses remains almost constant with Wi. Thus, the introduction of polymers increases the drag indirectly through increase in viscous stresses. The viscous drag in the x1 direction comes from contributions associated with the viscous 1 k stresses s11, s12 and s13 where s1k ¼ @u þ @u . The drag contributions @xk @x1 due to the viscous stresses are given by:

Cd/C0

drag due to

s11

drag due to

þ 4b s13 n3 dS S |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} drag due to

θ=0.315 θ=0.284

1.3

θ=0.351 θ=0.331

1.2

θ=0.384

1.1

1

Z Z Z Z ¼ 4b s11 n1 dS þ 4b s12 n2 dS S S |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} drag

C d;v isc |fflffl{zfflffl}

Z Z

1.4

0

0.5

1

1.5

Wi Fig. 7. Experimental Comparison. s, experimental [4]; j, simulation(comparison to [4]); —, Theory (h 1) [19]; }, simulation(comparison to [19]).

ð11Þ

polymer drag

viscous

1.5

viscous drag

s12

ð12Þ

s13

The contributions from these components with respect to the contributions at Wi = 0 are shown in Fig. 8d. The dominant contribution comes from the s11 component of the viscous stress. Thus, the presence of polymers in the fluid indirectly increases the drag on the sphere through the s11 component of the viscous stress. The drag coefficients of the simulations used to compare with van den Brule and Gheissary [4] and Housiadas and Tanner [19] are broken down into individual components using Eqs. (11) and (12). This is shown in Fig. 9. It is important to note that these results are presented for small Wi. The drag contributions due to all the three components increases with Wi. The increase in the drag contribution due to the viscous stresses is comparable to the increase due to the other two components for Wi 6 1.5. The results for Wi = 2.34 in Fig. 8c also suggests that there is an initial increase for all three components and the increase due to the viscous

55

S. Padhy et al. / Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60

1.5

3 2.5

1.3

Cd/C0

Cd/C0

1.4

1.2

θ=0.169

2

θ=0.227 1.5

1.1 1

1 0

2

4

6

8

10

5

10

Wi

(a)

(b)

15

Viscous Components

3

2 1.5 1 0.5

0

Wi

2.5

Drag Components

θ=0.123

θ=0.1

0

2

4

6

8

2.5 2 1.5 1 0.5

10

0

2

4

Wi

6

8

10

Wi

(c)

(d)

Fig. 8. (a) Effect of viscoelasticity on drag coefficients (h = 0.1). h, simulation. (b) Experimental comparison. s, experimental [5]; h, simulation (Sc = 106); M, simulation (Sc = 1). (c) Components of drag(h = 0.1). , polymer drag; , form drag; , viscous drag (d) Components of viscous drag(h = 0.1). - -, drag due to s12; – –, drag due to s11; - € -, drag due to s13.

1.4

1.4

Viscous Components

Drag Components

1.5

1.3 1.2 1.1 1 0.9

0

0.5

1

1.5

1.3 1.2 1.1 1 0.9

0

0.5

1

Wi

Wi

(a)

(b)

1.5

Fig. 9. Drag components for small Wi. h, simulation (comparison to [4]); s, simulation (comparison to [19]). (a) Components of drag. --, polymer drag;     , form drag; —-, viscous drag. (b) Components of viscous drag. --, drag due to s12; ——, drag due to s11; - - - -, drag due to s13.

stresses dominates for Wi P 2. The contributions from s11,s12 and s13 component of the viscous stress are of the same order for Wi 6 1.5. This is in accordance with the result for Wi = 2.34 in Fig. 8d as the contribution from the s11 component of the viscous stress becomes dominant for Wi > 2.34. Thus, for small Wi all the components play an equal role in increasing the drag on the sphere but as Wi increases the viscous stresses begin to dominate. In order to further understand the mechanism of interaction between the polymers and the sphere, the flow-induced stretch of the polymers and the s11 component of the viscous stress need to be examined. An isosurface of the trace of the conformation tensor for the polymers is plotted in Fig. 10. It can be seen that the polymers are stretched to large distances (of O(10D)) in the direction of the shear flow. Thus, any accurate simulation for a single sphere requires a large domain in the direction of shear flow for capturing these elongated regions. A domain size of 37.5D in the shear direction was used in our simulations because all the regions

where the trace of the conformation tensor is greater than the value at the inlet were contained inside the domain for Wi 6 15. The plots of the s11 component of the viscous stress on the sphere surface are shown in Fig. 11. Stress distributions for the case with only orthogonal shear flow (Fig. 11a and b) and the case with both uniform flow and orthogonal shear flow at Wi = 4.69 (Fig. 11c and d) are shown. We shall refer to the part of the sphere in the positive half-space in the x1 direction as the ‘‘front surface’’ of the sphere and the part in negative half-space as the ‘‘back surface’’. In the case with only orthogonal shear flow, there is positive stress on the top right and the bottom left on the front surface of the sphere and negative stresses on the rest of the front surface of the sphere. The stresses switch sign for the back surface of the sphere. The positive stresses are symmetric with respect to the negative stresses and thus they balance to give a net drag of zero. In the case with both uniform flow and orthogonal shear flow, a similar stress distribution is observed but now the stress

56

S. Padhy et al. / Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60

Fig. 10. Contour plots of trace of conformation tensor at Wi = 14.08 and h = 0.123. Isosurfaces at 50% of maximum value. In these plots, (x, y, z) corresponds to (x1, x2, x3), with x = x1 the direction of gravity and z = x3 the direction of the shear flow.

distribution is not symmetric which leads to a net drag in the x1 direction. Although the uniform flow is relatively weak (i.e. h = 0.1), this break in symmetry is enough to provide a net positive drag in the x1 direction which ultimately leads to the increase in

total drag that scales with Wi (i.e. depends on the stress distribution created by the orthogonal shear). The non-linear coupling of the shear and uniform motions accentuates this effect. Thus, a weak uniform flow combines non-linearly with the orthogonal shear flow by breaking the symmetry which ultimately leads to the increase in drag values. The contribution to the drag coefficients due to the polymer stresses is small compared to the viscous stresses. But this is still an additional component which is not present in the flows without the polymers. The dominant contribution comes from the sP13 component of the polymer stress. The stress distributions on the surface of the sphere for the sP13 component of the polymer stress is shown in the Fig. 12. Stress distributions for the case with only orthogonal shear flow (Fig. 12a and b) and the case with both uniform flow and orthogonal shear flow at Wi = 4.69 (Fig. 12c and d are shown. We shall refer to the part of the sphere in the positive half-space in the x3 direction as the ‘‘right side’’ of the sphere and the part in negative half-space as the ‘‘left side’’. The stress distributions are again symmetric for the case with only orthogonal shear flow. There is positive stress on the front right surface and negative stress on the back right surface and vice versa for the left surface.The introduction of the uniform flow breaks this symmetry which leads to a net positive drag in the x1 direction. The contribution due to the polymer stresses is small compared to the viscous stresses because the magnitude of the polymer stresses is smaller. Thus, the primary mechanism for the increase in the drag is the break in the symmetry of the

Fig. 11. Contour plots of s11 at Wi = 4.69. Only shear flow (a) front (b) back. Both uniform flow and shear flow (c) front (d) back.

S. Padhy et al. / Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60

Fig. 12. Contour plots of

57

sP13 at Wi = 4.69. Only shear flow (a) right side (b) left side. Both uniform flow and shear flow (c) right side (d) left side.

Fig. 13. Contour plots of s11 without polymers for combined uniform and shear flow (a) front (b) back.

polymer and viscous stress distributions on the sphere as the weak uniform flow is coupled with the orthogonal shear flow. The coupling between the two flows needs to be non-linear for the increase in the drag to be substantial. The non-linear coupling between the flows is only due to the presence of polymers. The stress distributions for the combined uniform and

shear flow without the presence of polymers is examined to test the effect of the non-linear coupling. The stress distributions on the sphere surface for the s11 component of the viscous stress are shown in Fig. 13. The stress distributions are much more symmetric than that which is observed in the case with polymers (Fig. 11c and d). In the absence of the polymers the weak

58

S. Padhy et al. / Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60

Fig. 14. Contour plot of 1 component of velocity (u1) on x1  x2 plane with x3 = 0.25D at Wi = 4.69 and h = 0.1. (a) with polymers and (b) without polymers.

Fig. 15. Contour plots of s11 on x1  x2 plane with x3 = 0. (a) Only uniform flow with polymers (Wi = 0 h = 0.1). (b) Combined uniform flow and shear flow with polymers (Wi = 4.69 h = 0.1). (c) Combined uniform flow and shear flow without polymers.

uniform flow combines linearly with the shear flow and therefore it is unable to break the symmetry. Thus, the increase in the drag is very small. The addition of the polymers leads to

the non-linear coupling of the two flows. This helps the weak uniform flow to break the symmetry and consequently increase the drag.

59

S. Padhy et al. / Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60

Viscous Components

3 2.5 2 1.5 1 0.5 0

2

4

6

8

10

Wi

(b)

(a)

(c) Fig. 16. (a) Spherical coordinate system. (b) Components of viscous drag in spherical coordinate system. (c) Contour plot of sr/ at Wi = 4.69.

It is interesting to examine the distribution of the u1 component of velocity since the s11 component of the viscous stress is the gradient of u1. The distribution of u1 for the cases with and without polymers is shown in Fig. 14. The velocity component, u1, varies over a wider range for the flow with polymers. The distributions are also different for the cases with and without polymers. Since the polymer stress term is coupled with the solution of flow field, the introduction of polymers modifies the flow field. This changes the distribution of s11 component of the viscous stress on the surface of the sphere which ultimately leads to the increase in drag coefficients. The distributions of s11 component of the viscous stress in the x1  x2 plane along with the streamlines are shown in Fig. 15. In the case with only the uniform flow in the presence of polymers (Fig. 15a), there are positive stresses downstream of the sphere and negative stresses upstream. As the orthogonal shearing motion is introduced (Fig. 15b), the streamline is stretched downstream of the sphere and relaxed upstream of the sphere. This leads to a very significant increase in the positive stresses (s11) downstream and decrease in the negative stresses (s11) upstream of the sphere. This change in stress distribution creates drag which oppose the motion of the sphere and ultimately reducing the settling rate. The distribution of s11 for the combined uniform and shear flow without the presence of polymers is shown in Fig. 15c. The streamlines are very different from the case where polymers are present as the introduction of polymers changes the flow field dramatically in particular bending the streamlines closer to the sphere. We now analyze the viscous and polymer stresses in a spherical coordinate system (shown in Fig. 16a). We can determine the role of shear and normal stresses using a spherical coordinate system because at each point the shear and normal stresses are aligned with the unit

, drag due to srw;

, drag due to sr/;

, drag due to srr.

normals in the spherical coordinate system. The components of viscous and polymer stresses are each transformed into the spherical coordinate system.The unit normals in r,w and / directions are represented by ^er ; ^ew and ^e/ respectively. The contribution to the drag coefficient due to the viscous stresses is broken into the individual contributions due to each of the stress components (srr, srw, sr/, sww, s// and sw/). The drag contribution due to the viscous stresses is given by:

Z Z ¼ 4b srr ^er sin w cos /dS S |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} drag

C d;v isc |fflffl{zfflffl}

viscous

drag due to

Z Z

srr

  þ 4b sr/ ^e/ sin w cos /  ^er sin / dS S |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} drag due tosr/

Z Z

  þ 4b srw ^ew sin w cos / þ ^er cos w cos / dS S |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} drag due to

Z Z

srw

þ 4b sww ^ew cos w cos /dS S |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Z Z

drag due to

sww

 4b s// ^e/ sin /dS S |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} drag due to

Z Z

s//

  þ 4b sw/ ^e/ cos w cos /  ^ew sin / dS S |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} drag due to

sw/

ð13Þ

60

S. Padhy et al. / Journal of Non-Newtonian Fluid Mechanics 197 (2013) 48–60

The contributions from these components with respect to the contributions at Wi = 0 are shown in Fig. 16b. Only the contributions from the major components are shown. The dominant contribution comes from the sr/ component of the viscous stress for Wi > 5. Thus, for large Wi, the viscous shear stress at each point on the sphere is the cause for the increase in the drag. The distribution of sr/ component of the viscous stress in the x1  x2 plane is shown in Fig. 16c. The viscous shear stress (sr/) is enhanced downstream of the sphere. The increase in the viscous shear stresses oppose the motion of the sphere by inducing asymmetric stress distributions on the surface of the sphere. In our simulations, we found that there is an increase in drag predicted for a single sphere which matches qualitatively with the experiments. This increase in drag is smaller than that measured in the experiments. It should be noted that these experiments were carried out in a relatively concentrated suspension (c = 0.2 by volume, Tonmukayakul et al. [5]). Earlier experiments by van den Brule and Gheissary [4] carried out for a single sphere had shown a smaller decrease in settling rate. As the simulations were conducted for a single sphere the increase in drag is smaller than that observed in the experiments with a relatively concentrated suspension. The effect of particle concentration is an excellent topic of future studies. 6. Conclusions Simulation of a sphere sedimenting in a cross shear flow of a FENE-P viscoelastic fluid has been performed. The shear flow is shown to reduce the settling rate, as seen experimentally. An increase in drag with increase in shear Weissenberg number is found, even though the shear viscosity decreases over this range of Wi, in agreement with the experiments. The increase in drag (or reduction in settling velocity) is due to viscous stresses, of which the s11 component is the dominant contribution. We recall that 1, 2 and 3 are, respectively, the directions of the uniform flow, the velocity gradient, and the shear velocity (u3 ¼ c_ x2 in this system); note that 1 would be the neutral or vorticity direction in the absence of sedimentation. Thus, the mechanistic basis for the increase in drag is that weak uniform flow due to the sedimentation breaks the symmetry seen for shear alone due to the non-linear coupling created by the presence of polymers. This allows the s11 component of the viscous stress to induce large changes to the drag by being asymmetric. The shearing motion changes the streamlines in such a way that the induced viscous stresses oppose the sedimentation of the sphere which leads to the decrease in settling rates. Acknowledgments The authors acknowledge the following award for providing computing resources that have contributed to the research results reported within this paper: MRI-R2: Acquisition of a Hybrid CPU/ GPU and Visualization Cluster for Multidisciplinary Studies in Transport Physics with Uncertainty Quantification (http://www. nsf.gov/awardsearch/showAward.do/?AwardNumber=0960306). This award is funded under the American Recovery and Reinvestment Act of 2009 (Public Law 111-5). In addition, this research is supported under a Stanford Graduate Fellowship. The authors SP, ESGS and GI also acknowledge financial support from Halliburton in order to complete the computational part of this study. The authors would like to thank Dr. D. Richter for help in debugging a part of the code used for simulation.

References [1] H.A. Barnes, A review of the rheology of filled viscoelastic systems, Rheol. Rev. 1 (2003) 1–36. [2] M.J. Economides, K.G. Nolte, Reservoir Stimulation, Prentice Hall, 1989. [3] J. Chatterji, J. Borchardt, Applications of water-soluble polymers in the oil field, J. Petrol. Technol. 33 (1981). [4] B.H.A.A. van den Brule, G. Gheissary, Effects of fluid elasticity on the static and dynamic settling of a spherical particle, J. Non-Newton. Fluid Mech. 49 (1993) 123–132. [5] N. Tonmukayakul, J.F. Morris, R. Prud’homme, Method for Estimating Proppant Transport and Suspend Ability of Viscoelastic Liquids, US Patent Application Publication Number US20110219856 (A1) (2011). [6] J. Mewis, N.J. Wagner, Current trends in suspension rheology, J. Non-Newton. Fluid Mech. 157 (2009) 147–150. [7] F. Snijkers, G.D. Avino, P.L. Maffetttone, F. Greco, M.A. Hulsen, J. Vermant, Rotation of a sphere in a viscoelastic liquid subjected to shear flow. Part II. Experimental results, J. Rheol. 53 (2009) 459–480. [8] G.D. Avino, M.A. Hulsen, F. Snijkers, J. Vermant, F. Greco, P.L. Maffetttone, Rotation of a sphere in a viscoelastic liquid subjected to shear flow. Part I. Simulation results, J. Rheol. 52 (2008) 1331–1346. [9] F. Snijkers, G.D. Avino, P.L. Maffetttone, F. Greco, M.A. Hulsen, J. Vermant, Effect of viscoelasticity on the rotation of a sphere in shear flow, J. Non-Newton. Fluid Mech. 166 (2011) 363–372. [10] K.D. Housiadas, R.I. Tanner, The angular velocity of a freely rotating sphere in a weakly viscoelastic matrix fluid, Phys. Fluids 23 (2011). 051702-1–4. [11] N. Goyal, J.J. Derksen, Direct simulations of spherical particles sedimenting in viscoelastic fluids, J. Non-Newton. Fluid Mech. 183–184 (2012) 1–13. [12] R. Scirocco, J. Vermant, J. Mewis, Effect of the viscoelasticity of the suspending fluid on structure formation in suspensions, J. Non-Newton. Fluid Mech. 117 (2004) 183–192. [13] D. Won, C. Kim, Alignment and aggregation of spherical particles in viscoelastic fluid under shear flow, J. Non-Newton. Fluid Mech. 117 (2004) 141–146. [14] M.K. Lyon, D.W. Mead, R.E. Elliot, L.G. Leal, Structure formation in moderately concentrated viscoelastic suspensions in simple shear flow, J. Rheol. 45 (2001). [15] R.J. Phillips, L. Talini, Chaining of weakly interacting particles suspended in viscoelastic fluid, J. Non-Newton. Fluid Mech. 147 (2007) 175–188. [16] S. Bobroff, R.J. Phillips, Nuclear magnetic resonance imaging investigation of concentrated suspensions in non-Newtonian fluids, J. Rheol. 42 (1998) 1419– 1436. [17] E. Allen, P.H.T. Uhlherr, Nonhomogenous sedimentation in viscoelastic fluids, J. Rheol. 33 (1989) 627–638. [18] S. Mora, L. Talini, C. Allain, Structuring sedimentation in a shear-thinning fluid, Phys. Rev. Lett. 95 (2005). 088301-1–4. [19] K.D. Housiadas, R.I. Tanner, The drag of a freely sendimentating sphere in a sheared weakly viscoelastic fluid, J. Non-Newton. Fluid Mech. 183–184 (2012) 52–56. [20] R.B. Bird, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, second ed., Wiley, 1987. [21] R. Sureshkumar, A.N. Beris, Effect of artificial stress diffusivity on the stability of numerical calculations and the flow dynamics of time-dependent viscoelastic flows, J. Non-Newton. Fluid Mech. 60 (1995) 53–80. [22] D. Richter, G. Iaccarino, E.S.G. Shaqfeh, Simulations of three-dimensional viscoelastic flows past a circular cylinder at moderate reynolds numbers, J. Fluid Mech. 651 (2010) 415–442. [23] H. Giesekus, A simple constitutive equation for polymer fluids based on the concept of deformation-dependent tensorial mobility, J. Non-Newton. Fluid Mech. 11 (1982) 69–109. [24] D. Richter, E.S.G. Shaqfeh, G. Iaccarino, Floquet stability analysis of viscoelastic flow over a cylinder, J. Non-Newton. Fluid Mech. 166 (2011) 554–565. [25] S. Kang, G. Iaccarino, F. Ham, Dns of buoyancy-dominated turbulent flows on a bluff body using the immersed boundary method, J. Comput. Phys. 228 (2009) 3189–3208. [26] K. Mahesh, G. Constantinescu, S. Apte, G. Iaccarino, F. Ham, P. Moin, Large-eddy simulation of reacting turbulent flows in complex geometries, J. Appl. Mech. 73 (2006) 374–381. [27] H. Yatou, Flow pattern transition accompanied with sudden growth of flow resistance in two-dimensional curvilinear viscoelastic flows, Phys. Rev. E 82 (2010) 036310(13). [28] D.O.A. Cruz, F.T. Pinho, P.J. Oliveira, Analytical solutions for fully developed laminar flow of some viscoelastic liquids with a newtonian solvent contribution, J. Non-Newton. Fluid Mech. 132 (2005) 28–35. [29] B. Yang, B. Khomami, Simulations of sedimentation of a sphere in a viscoelastic fluid using molecular based constitutive models, J. Non-Newton. Fluid Mech. 82 (1999) 429–452. [30] M.T. Arigo, D.R. Rajagopalan, N.T. Shapely, G.H. McKinley, The sedimentation of a sphere through an elastic fluid. Part 1. Steady motion, J. Non-Newton. Fluid Mech. 60 (1995) 225–257. [31] V. Levich, Physiochemical Hydrodynamics, Prentice Hall, 1962.