The Einstein viscosity with fluid elasticity

TL;DR

This paper develops a first-order correction to the viscosity of dilute suspensions in viscoelastic fluids, revealing shear-thickening behavior due to fluid elasticity and emphasizing the importance of extensional rheology modeling.

Contribution

It introduces a perturbation theory for suspension viscosity considering fluid elasticity effects up to second order in Weissenberg number, validated by numerical simulations.

Findings

Suspension shear-thickens near particles due to elastic stretching.

Extensional rheology significantly influences shear rheology predictions.

Strain-rate thickening observed at higher order in Weissenberg number.

Abstract

We give the first correction to the suspension viscosity due to fluid elasticity for a dilute suspension of spheres in a viscoelastic medium. Our perturbation theory is valid to $O(\phi\mathrm{Wi}^2)$ in the Weissenberg number $\mathrm{Wi}=\dot\gamma \lambda$, where $\dot\gamma$ is the typical magnitude of the suspension velocity gradient, and $\lambda$ is the relaxation time of the viscoelastic fluid. For shear flow we find that the suspension shear-thickens due to elastic stretching in strain hot spots near the particle, despite the fact that the stress inside the particles decreases relative to the Newtonian case. We thus argue that it is crucial to correctly model the extensional rheology of the suspending medium to predict the shear rheology of the suspension. For uniaxial extensional flow we correct existing results at $O(\phi\mathrm{Wi})$, and find dramatic strain-rate thickening at $O(\phi\mathrm{Wi}^2)$. We validate our theory with fully resolved numerical simulations.

Equations

$$\eta =(1+2.5\varphi +...)\mu ,$$

$$\overline{\mathbf{\sigma }(\mathbf{x})}=\int \mathbf{\sigma }(\mathbf{x}|\mathbf{y})P(\mathbf{y})\text{d}\mathbf{y}+O({\varphi }^2).$$

$$\overline{\mathbf{\sigma }}=\frac{\varphi }{V_p}{\int }_V\mathbf{\sigma }(\mathbf{r})\text{d}\mathbf{r}=\frac{\varphi }{V_p}{\int }_{V_p}\mathbf{\sigma }(\mathbf{r})\text{d}\mathbf{r}+\frac{\varphi }{V_p}{\int }_{V_f}\mathbf{\sigma }(\mathbf{r})\text{d}\mathbf{r},$$

$$\frac{\varphi }{V_p}{\int }_{V_p}\mathbf{\sigma }(\mathbf{r})\text{d}\mathbf{r}=\frac{\varphi }{V_p}{\int }_{S_p}\mathbf{r}(\mathbf{\sigma }\cdot \mathbf{n})\text{d}S=\frac{\varphi }{V_p}\mathbf{S}.$$

$$\nabla \cdot \mathbf{\sigma }=0,\nabla \cdot \mathbf{u}=0,$$

$$\mathbf{\sigma }=-p\mathbf{\delta }+2(1-{\mu }_r)\mathbf{e}+{\mu }_r\mathbf{\Pi },$$

$$\mathbf{\Pi }+\text{Wi}[(\mathbf{u}\cdot \nabla )\mathbf{\Pi }-\mathbf{a}\cdot \mathbf{\Pi }-\mathbf{\Pi }\cdot {\mathbf{a}}^{\mathrm{T}}]=\mathbf{a}+{\mathbf{a}}^{\mathrm{T}}.$$

$$\overline{\nabla \mathbf{u}}=\mathbf{E}+\mathbf{O}.$$

$$\overline{(\mathbf{u}\cdot \nabla )\mathbf{\Pi }}=0$$

$$\mathbf{u}\sim \overline{\mathbf{u}},\mathbf{\sigma }\sim \overline{\mathbf{\sigma }},|\mathbf{r}|\sim {\varphi }^{-1/3}.$$

$$\displaystyle\overline{\boldsymbol{\sigma}}=-\langle p\rangle_{F}\boldsymbol{\delta}\!+\!2\boldsymbol{E}+\frac{\phi}{V_{p}}\boldsymbol{S}+\mu_{r}\textrm{Wi}\big{\langle}\boldsymbol{a}\boldsymbol{\Pi}+\boldsymbol{\Pi}\boldsymbol{a}^{\rm T}\big{\rangle}_{F}\,,$$

$$⟨\mathbf{a}{⟩}_F\equiv \frac{\varphi }{V_p}{\int }_{V_f}\mathbf{a}\text{d}V.$$

$$\mathbf{u}={\mathbf{u}}^{(0)}+\text{Wi}{\mathbf{u}}^{(1)}+...,\mathbf{\Pi }={\mathbf{\Pi }}^{(0)}+\text{Wi}{\mathbf{\Pi }}^{(1)}+...,$$

$$p=p^{(0)}+\text{Wi}{p}^{(1)}+...,\mathbf{\omega }={\mathbf{\omega }}^{(0)}+\text{Wi}{\mathbf{\omega }}^{(1)}+....$$

$${\eta }_S=1+2.5\varphi +\left(0.62-0.03{\mu }_r\right)\varphi {\mu }_r{\text{Wi}}^2.$$

$${\alpha }_S^{\mathrm{stresslet}}$$

$${\alpha }_S^{\mathrm{fluid}}$$

$${\mathbf{S}}^A$$

$${\mathbf{S}}^B$$

$${\eta }_E=1+{\mu }_r\text{Wi}+3{\mu }_r{\text{Wi}}^2+2.5\varphi$$

$$+\varphi {\mu }_r\left(2.68\text{Wi}+9.36{\text{Wi}}^2-0.1{\mu }_r{\text{Wi}}^2\right)$$

$${\alpha }_E^{\mathrm{stresslet}}$$

$${\alpha }_E^{\mathrm{fluid}}$$

$${\eta }_S=1+\frac{5}{2}\varphi +\varphi {\mu }_r({\alpha }_S^{\mathrm{stresslet}}+{\alpha }_S^{\mathrm{fluid}})$$

$${\alpha }_S^{\mathrm{stresslet}}$$

$${\alpha }_S^{\mathrm{fluid}}$$

$${\eta }_E=1+{\mu }_r\text{Wi}+3{\mu }_r{\text{Wi}}^2+\frac{5}{2}\varphi +\varphi {\mu }_r({\alpha }_S^{\mathrm{stresslet}}+{\alpha }_S^{\mathrm{fluid}})$$

$${\alpha }_E^{\mathrm{stresslet}}$$

$${\alpha }_E^{\mathrm{fluid}}$$

$$\overline{\mathbf{a}}=\frac{\varphi }{V_p}{\int }_V\mathbf{a}(\mathbf{r})\text{d}V,$$

$$\overline{{\partial }_j{u}_i}=\frac{\varphi }{V_p}{\int }_V{\partial }_j{u}_i\text{d}V=\frac{\varphi }{V_p}{\int }_{S_{\infty }}{u}_i{n}_j\text{d}S.$$

$$\frac{\varphi }{V_p}{\int }_{S_{\infty }}(E_{ik}+O_{ik})r_k{n}_j\text{d}S=E_{ij}+O_{ij},$$

$${\int }_S{n}_i{r}_j\text{d}S={\int }_V{\partial }_i{r}_j\text{d}V=V{\delta }_{ij}$$

$$u_i$$

Full text

I Introduction

The rheology of complex suspensions is of fundamental importance in many fields of science and engineering. For most applications it not desirable, or practical, to resolve the details at the scale of the polymers or particles in the suspension. A fundamental problem in rheology is the development of methods to coarse-grain the microscopic fluid dynamics and thus create a useful macroscopic continuum description of the stress-strain response of a complex suspension. This is a difficult problem at the intersection of fluid dynamics and statistical mechanics.

A cornerstone of our theoretical understanding is so-called dilute suspension rheology Batchelor (1970). This approximation describes the stresses that arise due to isolated particles in a suspending medium, and particle-particle interactions are neglected. Mathematically, it arises at first order in a perturbation theory in the volume fraction $\phi$ of particles in the suspension Hinch (1977).

Einstein (1906, 1911) first devised an early form of dilute suspension rheology in his doctoral thesis. He showed that the bulk shear viscosity $\eta$ of a suspension of rigid, inertia-free, and neutrally buoyant spheres in a Newtonian fluid is

[TABLE]

where $\mu$ is the viscosity of the suspending fluid in the absence of any particles. The suspension viscosity increases because the particles resist deformation. Their internal stresses increase, which results in an $O(\phi)$ increase in the suspension viscosity. In a seminal paper, Batchelor (1970) described how to generalize Einstein’s calculation to compute the complete suspension stress tensor. He termed the increased particle stress the ‘stresslet’ contribution, because of its relation to the symmetric first moment of surface tractions over the particle surface.

When the suspending medium is viscoelastic, for example by the addition of polymers, there are two mechanisms that change the resulting suspension stress. First, the stresslet contribution may change, because the surface tractions change. Second, in contrast to the Newtonian case, there is additional stress in the fluid phase due to the polymers stretching in the flow gradients induced by the particle. We call this latter contribution the ‘particle induced fluid stress’ Yang et al. (2016).

Experiments with suspensions of spherical particles in viscoelastic fluids show shear-thickening at low ($<10\%$) volume fractions of particles Zarraga et al. (2001); Scirocco et al. (2005); Dai et al. (2014). In contrast, Newtonian suspensions shear-thin or thicken only at rather high volume fractions $\phi>30$-$40\%$, because of particle interactions Wagner and Brady (2009). These observations indicate that dilute suspension rheology may be useful to understand the mechanism of shear thickening of viscoelastic suspensions.

The first deviation from a Newtonian fluid for a viscoelastic fluid medium is described by the Second-order fluid Larson , which is an approximation in slow flows or, equivalently, fast relaxation of the elastic fluid. It is valid to linear order in the Weissenberg number $\textrm{Wi}=\dot{\gamma}\lambda$, where $\dot{\gamma}$ is the typical magnitude of the fluid velocity gradient, and $\lambda$ is the relaxation time of the viscoelastic fluid. Koch and Subramanian (2006, 2008) first calculated the correct rheology of a dilute suspension of spheres in a Second-order fluid, after several earlier attempts with conflicting results Kaloni and Stastna (1983); Sun and Jayaraman (1984); Mifflin (1985); Li and Slattery (1989); Greco et al. (2005, 2007a, 2007b); Housiadas and Tanner (2009). Recently, Yang et al. (2016) could discriminate between the different theoretical calculations by fully resolved numerical simulations. But in the Second-order fluid limit there is no correction to the shear viscosity. While there are Second-order fluid corrections to the normal stress differences due to fluid elasticity Koch and Subramanian (2006, 2008), Einstein’s result (1) remains the leading order correction to the suspension viscosity.

Yang et al. (2016) and Koch et al. (2016) independently studied suspension stress in shear flow by numerical simulation and a semi-analytical theory, respectively. They both found that as Wi is increased from vanishingly small values, the stresslet contribution to the dilute suspension shear viscosity decreases, but the particle induced fluid stress increases. Thus both found that the net effect is shear thickening of viscosity and the first normal stress.

In this paper we analytically calculate the suspension stress for any linear motion of the suspension, by a perturbation theory to $O(\phi\textrm{Wi}^{2})$. As a particular case we find the correction to the Einstein viscosity (1) due to fluid elasticity. Our calculation reveals how shear-thickening arises from strain ‘hot-spots’ in the disturbance flow around particles. In addition we analyze the stress in extensional flow, which is a fundamental rheological flow, and important in applications where, for example, a suspension is injected into a mold.

II Theory

II.1 Dilute suspension rheology

We consider an inertia-free suspension of rigid spherical particles of radius $a$ in a viscoelastic medium. We take the macroscopic flow to be a linear flow $\boldsymbol{U}=\boldsymbol{E}\cdot\boldsymbol{x}+\boldsymbol{O}\cdot\boldsymbol{x}$, where the rate of strain $\boldsymbol{E}$ is symmetric and the vorticity tensor $\boldsymbol{O}$ is anti-symmetric. The linear flow includes both the simple shear and extensional flows that are fundamental to rheology. In the following we first explain the general averaging procedure, we then give the details of the microscopic problem, and how to calculate the average.

The macroscopic stress at a point $\boldsymbol{x}$ in the suspension is an ensemble average of the microscopic stress $\boldsymbol{\sigma}(\boldsymbol{x})$, taken over all possible configurations of the suspension Batchelor (1970); Hinch (1977). In a dilute suspension it is unlikely that two particles are close enough to interact, and for low volume fraction $\phi$ the ensemble average is approximated by Hinch (1977)

[TABLE]

Here $\boldsymbol{\sigma}(\boldsymbol{x}\,|\,\boldsymbol{y})$ is the stress at $\boldsymbol{x}$ conditioned on the presence of a sphere centered at $\boldsymbol{y}$, and $P(\boldsymbol{y})$ is the probability to find a sphere at $\boldsymbol{y}$. Eq. (2) embodies the fact that we expect to find only a single particle within a volume $V\sim\phi^{-1}$. In a spatially homogenous suspension the probability $P(\boldsymbol{y})=\phi/V_{p}$ is uniform, where $V_{p}$ is the volume of a sphere, and the disturbance fields due to the presence of a particle depend only on $\boldsymbol{r}=\boldsymbol{x}-\boldsymbol{y}$.

It follows that the correct microscopic problem to solve, at this order, is that of a single particle centered at $\boldsymbol{r}=0$ in an asymptotically large volume of radius $R\sim\phi^{-1/3}$. Outside this volume we expect to find additional particles, and the assumptions that led to Eq. (2) are invalid. The problem formulation must be closed by a far field boundary condition for $|\boldsymbol{r}|\sim R$ that ensures a self-consistent theory, which we introduce after giving the equations of motion.

The ensemble average (2) thus reduces to

[TABLE]

where $V$ is the volume of the domain, including $V_{p}$, and the fluid volume $V_{f}=V-V_{p}$. The stress $\boldsymbol{\sigma}$ inside the particle is unknown, but via the divergence theorem and continuity of stress it is given by the stresslet of a freely suspended particle Batchelor (1970)

[TABLE]

II.2 Governing equations

The microscopic flow $\boldsymbol{u}$ is governed by the momentum equation and the incompressibility condition

[TABLE]

where $\boldsymbol{\sigma}$ is the stress tensor field. In the following all variables are non-dimensionalized: $\boldsymbol{r}^{\prime}=\boldsymbol{r}/a$, $t^{\prime}=\dot{\gamma}t$, $\boldsymbol{u}^{\prime}=\boldsymbol{u}/(\dot{\gamma}a)$, $\boldsymbol{\sigma}^{\prime}=\boldsymbol{\sigma}/(\mu\dot{\gamma})$, where $\mu$ is the shear viscosity at $\textrm{Wi}=0$ and $\dot{\gamma}=\sqrt{2\textrm{Tr}\,\boldsymbol{E}\boldsymbol{E}}$. We subsequently drop the primes from the notation. We use the Oldroyd-B constitutive model for the stress in a viscoelastic fluid. It represents a thermal bath of entropic springs, modeling polymers, that are transported and stretched by the fluid. The model captures the rheology of an elastic fluid without shear-thinning. The steady Oldroyd-B equations read Larson

[TABLE]

Here $p$ is the pressure, $\mu_{r}\boldsymbol{\Pi}$ is the stress due to the elastic polymers, and $\boldsymbol{a}$ is the flow gradient tensor with elements $a_{ij}=\partial u_{i}/\partial r_{j}$. The strain tensor is defined as $\boldsymbol{e}=(\boldsymbol{a}+\boldsymbol{a}^{\rm T})/2$, and the vorticity tensor $\boldsymbol{o}=(\boldsymbol{a}-\boldsymbol{a}^{\rm T})/2$, so that $\boldsymbol{a}=\boldsymbol{e}+\boldsymbol{o}$. The constitutive model has two parameters: the relaxation time $\lambda$ that appears in the Weissenberg number $\textrm{Wi}=\dot{\gamma}\lambda$, and the relative concentration of polymers given by the ratio $\mu_{r}=\mu_{p}/(\mu_{s}+\mu_{p})$ of the solvent ($\mu_{s}$) and polymer ($\mu_{p}$) contributions to the shear viscosity $\mu=\mu_{s}+\mu_{p}$ at $\textrm{Wi}=0$. The center-of-mass velocity $\boldsymbol{v}$ and the angular velocity $\boldsymbol{\omega}$ of the particle are determined by the condition that it is force- and torque-free.

II.3 Boundary conditions

The boundary condition on the surface of the sphere is the usual no-slip condition. The asymptotic far field boundary condition when $|\boldsymbol{r}|\sim R$ is determined by two self-consistency conditions. First, the suspension flow is linear by assumption, so the microscopic model must satisfy

[TABLE]

This implies that the flow field settles down to its mean value as $|\boldsymbol{r}|\sim\phi^{-1/3}$. Any correction due to the presence of another particle is of higher order in $\phi$ Hinch (1977). Second, the ensemble average

[TABLE]

in a homogenous suspension Koch et al. (2016); Rallison (2012) (see also Appendix B.) This implies that also the stress field settles down to its mean value as $|\boldsymbol{r}|\sim\phi^{-1/3}$. In summary,

[TABLE]

This argument is a mean field theory in the sense that we assume a macroscopic average field, and then determine the microscopic model to satisfy this assumption upon averaging. We demonstrate in Appendix B that these conditions indeed give a consistent theory.

Viscous flow problems in unbounded domains are often complicated by the slow algebraic decay of the disturbance fields, giving unphysical or divergent integrals. The self-consistency condition is in effect a regularization of those integrals, and our argument is akin to the regularization used to calculate the $O(\phi^{2})$ rheology in a Newtonian suspension O’Brien (1979). They considered an unbounded fluid domain, but instead invoked the asymptotic properties of $\boldsymbol{u}$ and $\boldsymbol{\sigma}$ when evaluating the integrals corresponding to our Eq. (3). We show in Appendix B that this procedure gives the same result as ours.

If one chooses to approximate the microscopic problem with an unbounded flow without regularization, it appears that the volume averaged stress is not equal to the ensemble averaged stress, even in a statistically homogenous suspension Koch and Subramanian (2006); Rallison (2012); Koch et al. (2016). This apparent ergodicity breaking has led to confusion regarding the correct averaging procedure, where some terms required ensemble averaging, whereas others could be volume averaged Greco et al. (2005); Koch and Subramanian (2006); Rallison (2012); Koch et al. (2016). We remove this ambiguity by correctly imposing the mean field conditions.

II.4 Averaging

It follows from Eq. (3), (4), and the Oldroyd-B constitutive equation (6) that

[TABLE]

where

[TABLE]

is short-hand for the integral over the fluid volume in Eq. (3). In Eq. (10) we used that the integrals $\langle\boldsymbol{e}\rangle_{F}=\boldsymbol{E}$, and $\langle(\boldsymbol{u}\cdot\nabla)\boldsymbol{\Pi}\rangle_{F}=0$ given the boundary conditions (9). In the following we omit any isotropic terms in the average stress, because they do not contribute to the suspension rheology. We denote the symmetric part of any tensor $\widehat{\boldsymbol{a}}\equiv(\boldsymbol{a}+\boldsymbol{a}^{\rm T})/2$. With the perturbation ansatz $\boldsymbol{\Pi}=\boldsymbol{\Pi}^{(0)}+\textrm{Wi}\,\boldsymbol{\Pi}^{(1)}+...$ we have

[TABLE]

where the stresslet $\boldsymbol{S}$ must be evaluated to $O(\textrm{Wi}^{2})$, the integrals in the first bracket must be evaluted to $O(\textrm{Wi})$, and those in the second bracket to $O(1)$.

The integrals over the quadratic and cubic terms are calculated by splitting $\boldsymbol{a}=\overline{\boldsymbol{a}}+\boldsymbol{a}^{\prime}$, and noting that $\langle\boldsymbol{a}^{\prime}\rangle_{F}=\phi\boldsymbol{E}$. This is because the average strain in the fluid phase is $(1+\phi)\boldsymbol{E}$, to compensate for the fact that there is no strain inside the particle. The remaining integrals of the type $\langle\boldsymbol{a}^{\prime}\cdot\boldsymbol{a}^{\prime}\rangle_{F}$, $\langle\boldsymbol{a}^{\prime}\cdot\boldsymbol{a}^{\prime}\cdot\boldsymbol{a}^{\prime}\rangle_{F}$, and $\langle\boldsymbol{a}^{\prime}\cdot[(\boldsymbol{u}\cdot\nabla)\boldsymbol{e}]\rangle_{F}$ are evaluated using the flow solutions that satisfy the boundary conditions (9) to $O(1)$ in $\phi$, which gives an accurate result to $O(\phi)$.

We evaluate the stresslet via the Lorentz reciprocal theorem Kim (1991); Einarsson and Mehlig (2017)

Here $\boldsymbol{\sigma}^{E}=\mu_{r}(\boldsymbol{\Pi}-2\boldsymbol{e})$ is the non-linear part of the stress tensor, and $\boldsymbol{M}$ is the rank three tensor such that $\boldsymbol{M}:\boldsymbol{E}$ is the Stokes solution for a sphere in a otherwise quiescent fluid, but with a strain flow $\boldsymbol{E}\cdot\boldsymbol{r}$ on the surface. See Appendix C for a detailed derivation of Eq. (LABEL:eqn:rt_stresslet).

We must know the flow and stress fields associated with the microscopic problem to $O(\textrm{Wi})$ to evaluate the average (LABEL:eqn:finalvolumeavg) and the stresslet (LABEL:eqn:rt_stresslet) to $O(\phi\textrm{Wi}^{2})$. The required flow and stress fields are given by a regular perturbation theory in Wi,

[TABLE]

We follow the method described in Ref. 25 to calculate these solutions. We perform the necessary algebra and integrations to evaluate Eqns. (LABEL:eqn:finalvolumeavg) and (LABEL:eqn:rt_stresslet) using the Matte software in Mathematica Jonas Einarsson (2017).

III Numerical methods

The numerical evaluations of the bulk stress is obtained by using a massively parallel flow solver based on an unstructured finite volume formulation. Details for the solver can be found in Yang et al. (2016) and the references therein. In those previous studies, the full Navier-Stokes equation was solved and the only modification we have made in this study is to remove the convective term such that the solver is an unsteady Stokes solver. At steady state, i.e. when the stresslet components and the particle-induced fluid stress components display a relative change of less than $0.001$ over a characteristic time scale, the solutions discretely satisfy the governing equations (5) and (6).

The computation domain is a cubic box with a sphere at the center. For the shear flow results, we use the same boundary conditions as reported in Ref. 5. At the sphere surface, the velocity corresponds to the solid body rotation of the sphere with angular velocity $\omega$, which is determined by iterating until the non-dimensionalized torque is less than 0.01. Equal and opposite velocities are applied on two walls of the computation box to drive the shear flow and periodic boundaries are applied on the remaining four walls. Note that we don’t need to apply boundary conditions for the polymer stress in this case because the two walls and sphere surface have zero mass flux. For the extensional flow results, we impose zero velocity on the sphere surface and $\boldsymbol{u}=\boldsymbol{E}\cdot\boldsymbol{r}$ on the computational box. In this case, since there is mass flux through the computation box, we also impose the ensemble averaged value of the polymer stress at the computation box to match the boundary conditions (9).

We have tested for mesh, time, and domain size convergence. Results change by less than 1.5% when we decrease the mesh size by a factor of 2, decrease the time step by a factor of 5, and increase the domain size by a factor of 2. We use a computation box that is 12 times the particle size. The mesh consists of tetrahedral elements that are finer on the sphere surface (0.03 particle diameters) and coarser on the boundary of the computation box (0.5 particle diameters) for a total of 2.5 million control volumes.

To extract the dependence of the stresslet on the value of the Weissenberg number we must subtract the $O(1)$ Newtonian contribution. Despite the numerical results for the viscosities being within $1.5$% of the theoretical Newtonian values at small values of Wi, subtracting the theoretical value is not precise enough to make a power law on a log-log scale. We therefore fit a curve $\eta=a+b\textrm{Wi}^{2}$ (or $a+b\textrm{Wi}+c\textrm{Wi}^{2}$ for extensional flow) to the three data-points with lowest Wi, and use $a$ as the value to subtract. The values of $a$ are within $1.5$% of the theoretical Newtonian value for all data-sets shown.

IV Results

We give our results in rounded decimal form for a concise and useful presentation. The exact expressions may be found in Appendix A.

Figure 1: Comparison of the theoretical results (lines) to fully resolved numerical simulations (markers) of the particle contribution to the suspension viscosities, as function of Wi. Top row shows suspension shear viscosity, bottom extensional viscosity in uniaxial strain. Left panels show stresslet contribution, right panels particle induced fluid stress. In all panels markers are numerical results for $\mu_{r}=0.68$ (red crosses), $\mu_{r}=0.5$ (green squares), and $\mu_{r}=0.01$ (blue circles). The solid lines show the $O(\textrm{Wi}^{2})$ theoretical results from Eq. (14) and Eq. (17). The dotted lines show the $O(\textrm{Wi})$ terms of Eq. (17).

IV.1 Shear viscosity

When the suspension flow is a simple shear flow $\boldsymbol{U}=y\hat{\boldsymbol{x}}$, we find that the shear viscosity $\eta_{S}=\overline{\sigma_{xy}}$ is

[TABLE]

The suspension is shear-thickening. The thickening is a consequence of two competing mechanisms. To show this we consider the components of the $O(\phi\textrm{Wi}^{2})$ contribution: $\eta_{S}=1+2.5\phi+\phi\mu_{r}(\alpha_{S}^{\mathrm{stresslet}}+\alpha_{S}^{\mathrm{fluid}})$.

[TABLE]

There is a negative, shear-thinning, contribution from the stresslet component, meaning that the stress in the sphere decreases due to the elasticity of the suspending fluid. The angular velocity also decreases at $O(\textrm{Wi}^{2})$ Housiadas and Tanner (2011): $\omega_{z}=-1/2+\mu_{r}\textrm{Wi}^{2}/4$. Thus both the symmetric and anti-symmetric moments of the surface tractions decrease at this order. On the other hand, the average polymer stress in the fluid phase gives a strong positive contribution to the shear viscosity. Our interpretation is that the polymers ‘absorb’ some of the stress that would otherwise have gone to surface traction on the particle. But that results in stretched polymers in the fluid around the particle that, on average, incur an even larger extra stress in the suspension. The net result is that the suspension shear-thickens. We compare our analytical result to fully resolved numerical simulations in Fig. IV (top row).

We further elucidate the mechanism of thickening by splitting the gradient $\boldsymbol{a}=\boldsymbol{e}+\boldsymbol{o}$ in Eq. (LABEL:eqn:finalvolumeavg). The resulting strain and vorticity correlations represent contributions from different flow types to the $O(\phi\textrm{Wi}^{2})$ viscosity correction. We present the contributions from all integrals, as well as the stresslet, in Fig. 2, for $\mu_{r}=1$. It is clear that the shear-thickening is due to the particle induced fluid stress in regions of flow around the particle with a strong straining component, represented by the integral $\langle\boldsymbol{e}\cdot\boldsymbol{e}\cdot\boldsymbol{e}\rangle_{F}$. The gradients in the regions fore and aft of the particle are strain-dominated, in the sense that the gradient tensor has three real eigenvalues. In the shear direction, however, the gradients are mixed strain and rotation. Regions of mixed flow diminish the effect of the strain, represented by the negative contributions of for example $\langle\boldsymbol{e}\cdot\boldsymbol{o}\cdot\boldsymbol{o}\rangle_{F}$. Nevertheless, the net result is thickening due to the regions of strain-dominated flow. We demonstrate this correlation in Fig. 3 which displays the local flow type side by side with the local contribution to the $O(\phi\textrm{Wi}^{2})$ particle induced fluid stress. The flow type is represented by the eigenvalue discriminant of the gradient tensor Chong et al. (1990), and the stress contribution is the integrand in the second bracket of Eq. (LABEL:eqn:finalvolumeavg).

We compare our theory to the semi-analytical results by Koch et al. (2016) for small values of Wi. We digitized the data from their Figs. 1 and 2 and compare to the relevant expressions from our theory in Fig. 4. The first two panels show the two contributions to the stresslet,

[TABLE]

The third panel shows the particle induced fluid stress. While the digitized data is too crude to make a quantitative comparison to a $\textrm{Wi}^{2}$ power law, we see that it is in good qualitative agreement with the theory.

The calculation by Koch et al. (2016) is valid to linear order in $\mu_{r}$. This approximation neglects the feedback mechanism that the elastic stresses modify the flow velocity field, which in turn generates new elastic stresses, which is represented by terms proportional to $\mu_{r}^{2}$ in our calculation. Our result shows that their model is a good approximation to $O(\phi\mu_{r}\textrm{Wi}^{2})$, because the terms of $O(\phi\mu_{r}^{2}\textrm{Wi}^{2})$ are very small.

Finally, we note that the normal stress differences are unaffected at $O(\phi\textrm{Wi}^{2})$, and at $O(\phi\textrm{Wi})$ we recover the known results Koch et al. (2016); Yang et al. (2016).

IV.2 Extensional viscosity

For uniaxial elongational flow $\boldsymbol{U}=x\hat{\boldsymbol{x}}-\frac{1}{2}y\hat{\boldsymbol{y}}-\frac{1}{2}z\hat{\boldsymbol{z}}$ we find that the elongational viscosity $\eta_{E}=(\overline{\sigma_{xx}}-\overline{\sigma_{yy}})/3$ is

[TABLE]

In this flow the leading correction to the viscosity is $O(\phi\textrm{Wi})$. Our result at this order differs from that of Greco et al. (2005) because they mistakenly evaluated the average over the fluid volume in Eq. (3) to $O(1)$ instead of $O(\phi)$ Greco et al. (2007b); Rallison (2012). They therefore report only the stresslet contribution, which is in agreement with our result for the stresslet. We separate the contributions from the stresslet and the particle induced fluid stress as $\eta_{E}=1+\mu_{r}\textrm{Wi}+3\mu_{r}\textrm{Wi}^{2}+2.5\phi+\phi\mu_{r}(\alpha_{S}^{\mathrm{stresslet}}+\alpha_{S}^{\mathrm{fluid}})$:

[TABLE]

Both stresslet and particle induced fluid stress thicken the suspension, and the dominant contribution arises from the fluid stress around the particle. We note that the coefficient of $\textrm{Wi}^{2}$ is so large that the second order result is important already for $\textrm{Wi}\approx 0.1$ (Fig. IV, bottom row.)

V Conclusions

We have calculated the stress in a dilute suspension of rigid spherical particles in an elastic fluid to $O(\phi\textrm{Wi}^{2})$. This revealed the first effects of elasticity on the suspension viscosity. When shearing the suspension we found that the stress in the particle phase decreases relative to the Newtonian case, but in return the stress in the elastic fluid surrounding the spheres increases. The net result is that the suspension is shear thickening. The main contribution to this thickening is the enhanced polymer stress in the straining regions of the fluid surrounding the spheres. This is significant because it implies that it is crucially important to correctly model the extensional rheology of the suspending fluid, in order to model the shear rheology of the suspension.

The available measurements of shear rheology in viscoelastic particle suspensions display shear thickening Zarraga et al. (2001); Scirocco et al. (2005); Dai et al. (2014). It is, however, not possible to extract an $\textrm{Wi}^{2}$ trend from their data for quantitative comparison. Further, the suspending fluids in the experiments were characterized by their shear rheology alone. We now believe that a successful quantitative comparison to experiment requires modeling the extensional rheology of the suspending medium. Thus, we argue that experimental measurements correlating the suspension rheology to the extensional rheology of the suspending medium will be valuable to further our understanding of complex suspensions at all values of Wi.

The contributions proportional to $\mu_{r}^{2}\textrm{Wi}^{2}$ represent the elastic stress disturbing the flow field, which in turn generates new elastic stresses. This effect is very small at this order of approximation, and therefore the $O(\mu_{r})$ semi-analytical theory for low polymer concentration Koch et al. (2016) is a good approximation at least up to $O(\phi\textrm{Wi}^{2})$.

We also gave the first correction to the extensional rheology of the suspension. In this case both the stresslet and particle induced fluid stresses contribute to significant thickening of the suspension with increasing strain rate. A recent experiment Dai and Tanner (2017) shows strain-rate thickening in a rather dense suspension ($\phi=40\%$). We are, however, not aware of any experimental measurements of the extensional viscosity of a dilute suspension of spheres in an elastic fluid.

Appendix A Exact results

In the main text we present our results in rounded decimal form for convenience. For reference, the exact results are for shear flow:

[TABLE]

[TABLE]

For uniaxial extensional flow:

[TABLE]

with

[TABLE]

Appendix B On the boundary conditions

Here we give some detail to the averaging integrals mentioned in the Letter. For the case of the average velocity gradient, and for the convective term in the polymer stress, we demonstrate that our calculation and the regularization scheme of O’Brien (1979) are equivalent. We also demonstrate explicitly that considering an unbounded domain without employing any regularization gives an inconsistent result.

In the following the full domain is denoted by $V$, and $S_{\infty}$ is a far field surface at the asymptotically large distance $R$ (taken to $\phi^{-1/3}$ in our Letter, but to $\infty$ in case of an unbounded domain.) The ensemble average of a quantity $\boldsymbol{a}$ is, in the dilute approximation,

[TABLE]

of which Eqns. (3) and (4) in the letter display the case of the stress tensor.

B.1 Average of the velocity gradient

The average $\overline{\partial_{j}u_{i}}=E_{ij}+O_{ij}$ by assumption in the mean field theory. The direct calculation is, by the divergence theorem,

[TABLE]

Our boundary condition requires $u_{i}\sim(E_{ik}+O_{ik})r_{k}$, when $r\sim\phi^{-1/3}$:

[TABLE]

where we used that $\phi/V_{p}=1/V$, and

[TABLE]

for any closed domain $V$ bounded by the surface $S$.

The method of O’Brien (1979) amounts to arguing that in the differential surface element $\textrm{d}S$, $\boldsymbol{r}$ samples a large enough region to achieve statistical stationarity for $\boldsymbol{u}$, and therefore $\boldsymbol{u}$ may be replaced by its statistical average in the integral. This obviously yields the same integral as our boundary condition.

Finally, should one assume an unbounded domain approximation for $\boldsymbol{u}$, the leading order flow field is

[TABLE]

one finds

[TABLE]

Here we took $S_{\infty}$ to be a sphere at distance $R$, as $R\to\infty$. The integral is convergent as $R\to\infty$ because $\boldsymbol{u}$ decays at the same rate as the surface element $\textrm{d}S$ grows. However, the result does not satisfy the mean field theory to $O(\phi)$.

B.2 Average of the convective term

The general argument for why $\overline{(\boldsymbol{u}\cdot\nabla)\boldsymbol{\Pi}}=0$ is as follows. The ensemble average commutes with the gradient operator, and $\boldsymbol{u}$ is divergence-free, so

[TABLE]

Split $\boldsymbol{u}=\overline{\boldsymbol{u}}+\boldsymbol{u}^{\prime}$ and $\boldsymbol{\Pi}=\overline{\boldsymbol{\Pi}}+\boldsymbol{\Pi}^{\prime}$, where primed quantities denote fluctuations around the mean. Then

[TABLE]

Because $\boldsymbol{\Pi}$ and $\boldsymbol{u}^{\prime}$ are statistically homogenous in the spatial variables by assumption, the gradient of their averages vanish. This argument was given earlier by Koch and Subramanian (2006) and Rallison (2012).

The direct calculation from the microscopic problem proceeds by the divergence theorem:

[TABLE]

The polymer stress vanishes inside the particle, and $u_{k}n_{k}=0$ on the particle surface. Our boundary condition requires $u_{k}\sim(E_{kl}+O_{kl})r_{l}$, and $\Pi_{ij}\sim\overline{\Pi_{ij}}$, when $r\sim\phi^{-1/3}$:

[TABLE]

The tensor $\boldsymbol{\Pi}$ is spatially homogenous by assumption, so $\overline{\Pi_{ij}}$ is a constant. The integral vanishes because $E_{kk}=O_{kk}=0$.

The argument of O’Brien (1979) leads to the same conclusion. In that case we would argue that both $\boldsymbol{u}$ and $\boldsymbol{\Pi}$ are statistically stationary within the surface element $\textrm{d}S$, and may therefore be replaced by their ensemble averages in the integral.

Finally, let us evaluate the volume averaged convective term using the leading order flow solutions in an unbounded domain without any regularization. The far field asymptotes of those fields read

[TABLE]

After a bit of algebra it follows that

[TABLE]

The integrand decays as $1/r^{2}$, which gives a convergent, but inconsistent result for the integral. We note that the volume average may also be evaluated directly as a volume integral. In that case a scaling analysis indicates that the integral is divergent, because the integrand decays as $1/r^{3}$. But that contribution vanishes identically upon integration of the angular variables, and the remaining radial integrals in fact decay as $1/r^{6}$ or faster, yielding the result in Eq. (33).

Appendix C Reciprocal theorem

Let $\boldsymbol{u}$ and $\boldsymbol{\sigma}^{N}$ denote the flow field and Newtonian stress tensor of a flow that is governed by the equation of motion

[TABLE]

The general form of Lorentz reciprocal theorem is Kim (1991)

[TABLE]

where $\tilde{\boldsymbol{u}}$ and $\tilde{\boldsymbol{\sigma}}$ are the flow and stress fields of an auxiliary Newtonian flow problem in the same geometry but with different boundary conditions. The volume integral is over a volume of fluid, and the surface integrals are to be taken over all boundaries of $V$ with surface normals pointing out from $V$.

Take the auxiliary problem to be the Stokes flow governed by

[TABLE]

Then

[TABLE]

We apply this theorem to our problem by noting that

[TABLE]

where

[TABLE]

Thus

[TABLE]

The surface integral on the left hand side is related to the stresslet because $\tilde{u}_{i}=\tilde{E}_{ik}r_{k}$ on the surface of the particle. The volume integral is well defined because $\tilde{u}_{i}\to 0\,,~{}r\to\infty$. The auxiliary flow is thus the disturbance flow around a sphere in a straining flow, and it decays as $1/r^{2}$ as $r\to\infty$. It is given by

[TABLE]

Split the surface integrals into the integral over the particle surface $S_{p}$ and the far-field surface $S_{\infty}$ that we take to be a sphere at $r=R_{\infty}$. Using the boundary conditions for both $\boldsymbol{u}$ and $\tilde{\boldsymbol{u}}$ we have

[TABLE]

Here $\boldsymbol{\omega}$ is the angular velocity of the particle. Starting from the left have the sought stresslet

[TABLE]

Next,

[TABLE]

Here we used that $\boldsymbol{\sigma}^{N}\sim 2\boldsymbol{E}+O(1/r^{2})$ and $\tilde{\boldsymbol{u}}\sim\frac{5}{2}\boldsymbol{r}(\boldsymbol{r}\boldsymbol{r}:\tilde{\boldsymbol{E}})/r^{5}+O(1/r^{4})$, as $r\to\infty$. Next,

[TABLE]

is the stresslet due to the non-linear polymer surface traction. Next,

[TABLE]

where we used that

[TABLE]

Next,

[TABLE]

because the particle in the auxiliary problem is torque-free.

Taken together we have, as $R_{\infty}\to\infty$,

[TABLE]

Finally, by renaming the dummy indices and expressing $\tilde{u}_{l}=M_{lik}\tilde{E}_{ik}$ in the volume integral, we eliminate $\tilde{\boldsymbol{E}}$ and arrive at the theorem stated in the main manuscript:

[TABLE]

Appendix D Derivation of non-linear stress

In this Appendix we employ a shorthand for contraction of rank two tensors, so that $\boldsymbol{a}\cdot\boldsymbol{\Pi}$ is written simply $\boldsymbol{a}\boldsymbol{\Pi}$, in the interest of keeping the expressions brief. There are no outer products of rank two tensors in this manuscript.

Starting from the constitutive equation

[TABLE]

insert $\boldsymbol{\Pi}=\boldsymbol{\Pi}^{(0)}+\textrm{Wi}\,\boldsymbol{\Pi}^{(1)}+\textrm{Wi}^{2}\,\boldsymbol{\Pi}^{(2)}$ and compare order by order:

[TABLE]

Inserting Eq. (52) in Eq. (53), and that result into Eq. (54), find

[TABLE]

[TABLE]

Thus

[TABLE]

To be consistent to $O(\textrm{Wi}^{2})$ the flow fields in the first bracket must be evaluated to $O(\textrm{Wi})$, whereas the fields in the second bracket only require the $O(1)$ Stokes flow solution.