Fluctuating stresslets and the viscosity of colloidal suspensions

TL;DR

This paper reveals a previously unobserved discrepancy in colloidal suspension viscosity measurements, which is resolved by incorporating fluctuating Brownian stresses and non-fluctuating stresslets into the hydrodynamic model, ensuring consistency between different calculation methods.

Contribution

It introduces an extended fluctuation-dissipation theorem that accounts for fluctuating Brownian stresses and non-fluctuating stresslets, improving viscosity measurement consistency in colloidal suspensions.

Findings

Discrepancy between shear stress and Green-Kubo viscosity methods identified.

Inclusion of fluctuating Brownian stresses resolves the discrepancy.

Non-fluctuating stresslets are essential for consistent viscosity calculations.

Abstract

Theory and simulation of Brownian colloids suspended in an implicit solvent, with the hydrodynamics of the fluid accounted for by effective interactions between the colloids, are shown to yield a marked and hitherto unobserved discrepancy between the viscosity calculated from the average shear stress under an imposed shear rate in the Stokesian regime and the viscosity extracted by the Green-Kubo formalism from the auto-correlations of thermal stress fluctuations in quiescent equilibrium. We show that agreement between both methods is recovered by accounting for the fluctuating Brownian stresses on the colloids, complementing and related to the traditional fluctuating Brownian forces and torques through an extended fluctuation-dissipation theorem based on the hydrodynamic grand resistance matrix. Time-averaging of the fluctuating terms gives rise to novel non-fluctuating stresslets. Brownian Dynamics simulations of spheroidal particles illustrate the necessity of these fluctuating and non-fluctuating contributions to obtaining consistent viscosities.

Equations

$${\eta }_s={\eta }_0\left(1+B\varphi \right),$$

$$\left(\begin{array}[]{c}\mathbf{F}_{\mathrm{H}}\\ \mathbf{S}_{\mathrm{H}}\end{array}\right)=-\left(\begin{array}[]{cc}\bm{\xi}^{F}_{U}&\bm{\xi}^{F}_{E}\\ \bm{\xi}^{S}_{U}&\bm{\xi}^{S}_{E}\end{array}\right)\left(\begin{array}[]{r@{\;-\;}l}\mathbf{U}&\mathbf{u}_{\infty}\\ &\mathbf{E}_{\infty}\end{array}\right),$$

$$\left(\begin{array}[]{c}\mathbf{U}\\ \mathbf{S}\end{array}\right)=\left(\begin{array}[]{cc}\bm{\mu}^{U}_{F}&\bm{\mu}^{U}_{E}\\ \bm{\mu}^{S}_{F}&\bm{\mu}^{S}_{E}\end{array}\right)\left(\begin{array}[]{c}\mathbf{F}_{\Phi}+\delta\mathbf{F}\\ -\mathbf{E}_{\infty}\end{array}\right)+\left(\begin{array}[]{c}\mathbf{u}_{\infty}\\ \mathbf{0}\end{array}\right),$$

$$⟨\delta \mathbf{F}(t)\otimes \delta \mathbf{F}(t^{\prime })⟩=2k_{B}T{\mathbf{\xi }}_U^F\delta (t-t^{\prime }),$$

$$B=\frac{1}{10k_{B}T{\eta }_{0}v}{\int }_0^{\infty }⟨\Delta \mathbf{S}(t):\Delta \mathbf{S}(0)⟩dt,$$

$$\left\langle\left(\begin{array}[]{c}\delta\mathbf{F}(t)\\ \delta\mathbf{S}(t)\end{array}\right)\otimes\left(\begin{array}[]{c}\delta\mathbf{F}(t^{\prime})\\ \delta\mathbf{S}(t^{\prime})\end{array}\right)\right\rangle=2k_{B}T\left(\begin{array}[]{cc}\bm{\xi}^{F}_{U}&\bm{\xi}^{F}_{E}\\ \bm{\xi}^{S}_{U}&\bm{\xi}^{S}_{E}\end{array}\right)\delta(t-t^{\prime}).$$

$$\left(\begin{array}[]{c}\mathbf{U}\\ \mathbf{S}\end{array}\right)=\left(\begin{array}[]{cc}\bm{\mu}^{U}_{F}&\bm{\mu}^{U}_{E}\\ \bm{\mu}^{S}_{F}&\bm{\mu}^{S}_{E}\end{array}\right)\left(\begin{array}[]{c}\mathbf{F}_{\Phi}+\delta\mathbf{F}\\ -\mathbf{E}_{\infty}\end{array}\right)+\left(\begin{array}[]{c}\mathbf{u}_{\infty}\\ \delta\mathbf{S}\end{array}\right).$$

$$\displaystyle\begin{split}\Delta\mathbf{Q}(t)&=\mathbf{Q}(t+\Delta t)-\mathbf{Q}(t)\\ &=-\bm{\mu}^{U}_{F}\bm{\nabla}_{\!Q}\Big{(}\Phi-\ln g_{Q}^{1/2}\Big{)}\Delta t-\bm{\mu}^{U}_{E}\mathbf{E}_{\infty}\Delta t\\ &\phantom{==}+\mathbf{u}_{\infty}\Delta t+k_{B}T\bm{\nabla}_{\!Q}\cdot\bm{\mu}^{U}_{F}\Delta t\\ &\phantom{==}+\bm{\mu}^{U}_{F}\left[\bm{\rho}^{F}_{U}\bm{\Theta}^{U}(t)+\bm{\rho}^{F}_{E}\bm{\Theta}^{E}(t)\right]\sqrt{\Delta t},\end{split}$$

$$\displaystyle\begin{split}\bar{\mathbf{S}}(t)&=\frac{1}{\Delta t}\int_{t}^{t+\Delta t}\mathbf{S}(t^{\prime})dt^{\prime}\\ &=\bm{\mu}^{S}_{F}\mathbf{F}_{\Phi}-\bm{\mu}^{S}_{E}\mathbf{E}_{\infty}+\bm{\nabla}_{\!X}\cdot\bm{\mu}^{S}_{F}\\ &\phantom{==}+\left[\left(\bm{\nabla}_{\!X}\bm{\rho}^{S}_{U}\right)\bm{\rho}^{F}_{U}+\left(\bm{\nabla}_{\!X}\bm{\rho}^{S}_{E}\right)\bm{\rho}^{F}_{E}\right]:\bm{\mu}^{U}_{F}\\ &\phantom{==}+\frac{1}{\sqrt{\Delta t}}\Big{\{}\bm{\mu}^{S}_{F}\left[\bm{\rho}^{F}_{U}\bm{\Theta}^{U}(t)+\bm{\rho}^{F}_{E}\bm{\Theta}^{E}(t)\right]\\ &\phantom{======}+\bm{\rho}^{S}_{U}\bm{\Theta}^{U}(t)+\bm{\rho}^{S}_{E}\bm{\Theta}^{E}(t)\Big{\}}.\end{split}$$

Full text

I Introduction

Einstein derived in his thesis that adding rigid spherical colloids to a Newtonian fluid of viscosity $\eta_{0}$ creates a suspension of effective viscosity

[TABLE]

with Einstein coefficient $B=5/2$, for low colloidal volume fractions $\phi$ Einstein1906 ; Einstein1911 . This celebrated result, based on the analytic solution of Stokesian straining flow around a spherical particle KimKarrila ; GuazzelliMorris , is readily reproduced by Brownian Dynamics (BD) simulations of an isolated spherical colloid suspended in a fluid subject to a linear shear flow. Viscosities can also be determined from quiescent fluids, using the Green-Kubo formalism of integrating the auto-correlations of the spontaneous stress fluctuations McQuarrie . Rather surprisingly, both the aforementioned theory and BD simulations of an isolated spherical particle then yield $B=0$, as will be demonstrated below. One might argue that this difference is an artefact of studying a one-particle system, which could be the reason that it appears not to have been discussed before in the literature, but we are of the opinion that it reveals a deficiency in the current understanding of stress calculations of suspensions of Brownian particles. We propose a solution, the inclusion of fluctuating Brownian stresses, that recovers agreement between equilibrium and non-equilibrium evaluations of the Einstein coefficient of a spherical particle, at $B=5/2$. By expressing the stochastic equations for the motion and stress in the Itô form, two novel non-vanishing stress contributions emerge from correlations between the various fluctuations. These terms affect the Einstein coefficients of isolated non-spherical particles, both in quiescent fluids and in flowing fluids. The novel terms also contribute to the viscosities of non-dilute solutions, including suspensions of spherical particles Batchelor72 ; Batchelor77 .

II Brownian motion

Consider a non-Brownian particle with generalized velocity $\mathbf{U}=(\mathbf{v},\bm{\omega})$, where $\mathbf{v}$ and $\bm{\omega}$ denote linear and angular velocity, respectively, in an incompressible fluid subject to a linear flow field characterized by the local generalized velocity $\mathbf{u}_{\infty}$ and the strain rate tensor $\mathbf{E}_{\infty}$, i.e. the traceless symmetric $(3\times 3)$ velocity gradient matrix. In the limit of Stokesian flow, the generalized hydrodynamic force $\mathbf{F}_{\mathrm{H}}$, a vector comprizing three force components and three torque components, and the deviatoric hydrodynamic stress $\mathbf{S}_{\mathrm{H}}$, a traceless symmetric $(3\times 3)$ matrix, acting on the particle are given by Durlofsky87 ; KimKarrila ; GuazzelliMorris

[TABLE]

where $\bm{\xi}$ is the grand resistance matrix KimKarrila ; GuazzelliMorris , whose four parts are labeled with a lower index specifying the multiplication partner and an upper index highlighting the ensuing result. Upon neglecting inertial effects, the equation of motion for a colloid experiencing also a generalized potential-derived force $\mathbf{F}_{\Phi}$ and a fluctuating Brownian force $\delta\mathbf{F}$ is solved from a balance of forces, $\mathbf{F}_{\Phi}+\mathbf{F}_{\mathrm{H}}+\delta\mathbf{F}=\mathbf{0}$. A partial inversion of the above equation then gives Durlofsky87 ; KimKarrila ; GuazzelliMorris

[TABLE]

where $\bm{\mu}$ is the generalized mobility matrix, again expressed as a combination of four labeled parts, and $\mathbf{S}=-\mathbf{S}_{\mathrm{H}}$ denotes the stress exerted on the fluid by the moving colloid. It is important to realize that the stress $\mathbf{S}$ is not the result of a balance of stresses, but a direct consequence of the velocity difference between colloid and fluid. Should one desire so, for instance when studying easily deformable particles, a balance must be constructed between the total stress acting on the colloid – due to the hydrodynamic stress $\mathbf{S}_{\mathrm{H}}$, a potential-derived $\mathbf{S}_{\Phi}$ and Brownian contributions (as discussed below) – and the internal elastic stress of the particle $\mathbf{S}_{\mathrm{I}}$; after solving this balance for the unknown $\mathbf{S}_{\mathrm{I}}$, the combination thereof with the elasticity tensor of the colloid yields its deformation. We will here consider rigid particles instead, and note for completeness that their stress balances are closed by unspecified Lagrange multipliers for $\mathbf{S}_{\mathrm{I}}$ at vanishing deformation. The random force perturbations $\delta\mathbf{F}$ have zero mean, are uncorrelated in time (Markovian) and obey the classical fluctuation-dissipation theorem derived from the symmetric positive-definite $(6\times 6)$ force-velocity segment of the resistance matrix Ermak78 ; Brady93 ,

[TABLE]

where $t$ and $t^{\prime}$ denote times, $k_{B}$ Boltzmann’s constant, $T$ the temperature and $\delta$ the Dirac delta function. The textbook proof of the fluctuation-dissipation theorem is its ability, in combination with a second-order Langevin equation of motion, to reproduce the Maxwell-Boltzmann equilibrium velocity distribution McQuarrie ; Kampen . Following the introduction of fluctuating hydrodynamics by Landau and Lifshitz Landau57 ; LifshitzPitaevskii , several authors have shown that the above theorem for a colloid also follows from the fluctuation-dissipation theorem of the fluid Zwanzig64 ; Fox70 ; Hauge73 ; Bedeaux74 ; Noetinger90 ; Singh17 . Note that the force perturbations affect the velocity difference between colloid and fluid and thereby give rise to an indirect Brownian contribution to the stress Bossis89 , $\delta\mathbf{S}_{\delta\mathbf{F}}=\bm{\mu}^{S}_{F}\delta\mathbf{F}(t)$, as follows from Eq. (3).

For a free spherical particle of volume $v$ in a linear shear flow, the above expressions give rise to the average stress $\langle\mathbf{S}\rangle=-\bm{\mu}^{S}_{E}\mathbf{E}_{\infty}$, where the minus sign indicates resistance to the flow. Inserting the theoretical expression for $\bm{\mu}^{S}_{E}$ KimKarrila ; GuazzelliMorris then yields $B=5/2$, as expected. But applying the Green-Kubo formalism to the spontaneous stress fluctuations $\Delta\mathbf{S}(t)=\mathbf{S}(t)-\langle\mathbf{S}\rangle$ in a quiescent fluid McQuarrie ; Daivis94 ,

[TABLE]

yields $B=0$ for a spherical particle, as follows from observing that under these conditions $\Delta\mathbf{S}(t)=\delta\mathbf{S}_{\delta\mathbf{F}}(t)=\mathbf{0}$ since $\bm{\mu}^{S}_{F}=\mathbf{0}$ for a spherical particle. It is this disconcerting discrepancy between the $B$-s that motivates the current research.

III Fluctuating Brownian stress

A suspended colloidal particle experiences a myriad of collisions due to the thermal motions of the surrounding solvent molecules. The sum over all collisions over a short time interval – sufficiently short to ignore the motion of the colloid, yet encompassing a large number of molecular collisions – gives rise to the fluctuating Brownian force $\delta\mathbf{F}(t)$ discussed above. Note that this generalized force comprises a force $\delta\mathbf{f}(t)$ and a torque $\delta\bm{\tau}(t)$, which represent distinct ‘projections’ of the same molecular noise integrated over the surface of the colloid Batchelor70 ; Durlofsky87 ; GuazzelliMorris . Their common origin implies that the force and torque are correlated, as reflected by the fluctuation-dissipation theorem in Eq. (4). It is only natural to assume that the collisions also give rise to a fluctuating Brownian stress on the particle, $\delta\mathbf{S}$, i.e. a stress distribution over the surface that would cause a soft particle to deform, which constitutes a third projection of the same molecular noise. This direct fluctuating Brownian stress $\delta\mathbf{S}$ is not to be confused with $\delta\mathbf{S}_{\delta\mathbf{F}}$, the latter being an indirect fluctuating stress resulting from a velocity difference between colloid and fluid caused by the first and second projections of the thermal noise. We next need to determine the strength of the fluctuating Brownian stress $\delta\mathbf{S}$, which in view of the preceding discussion does not follow from a stress balance on the particle. To conform with common practice in the field KimKarrila ; GuazzelliMorris , our interest here will be on the deviatoric parts of the stresses.

For any (non-Brownian) colloid experiencing a flow field, the hydrodynamic force acting on the particle is obtained as the zeroth moment of the traction vector integral over the surface while the torque and the stress or ‘stresslet’ on the particle are given by (a permutation of) the anti-symmetric and the symmetric first moments of the traction vector integral over the surface, respectively Batchelor70 ; Durlofsky87 ; GuazzelliMorris . If this flow field is replaced by the fluctuating hydrodynamics of the fluid, the strengths of the resulting fluctuating Brownian force $\delta\mathbf{f}$ and torque $\delta\bm{\tau}$, as well as their cross-correlations, are given by the fluctuation dissipation theorem of Eq. (4). We now hypothesize that the fluctuations of $\delta\mathbf{F}$ and $\delta\mathbf{S}$, given their common origin as projections of fluctuating hydrodynamics, are related by a generalized fluctuation-dissipation theorem based on the grand resistance matrix,

[TABLE]

Whereas the force correlations in this expression are well established Zwanzig64 ; Fox70 ; Hauge73 ; Bedeaux74 ; Noetinger90 ; Singh17 , the stress correlations and the force-stress cross-correlations have attracted little attention. A recent study supports the validity of our assumption for spherical particles, however without providing an explicit expression Singh17 . Introducing the fluctuating stresslet into Eq. (3) gives the extended expression for the motion and the stress,

[TABLE]

It is evident that this expression produces identical translational and rotational Brownian motion to the classical expression, even though the fluctuating force is now correlated to the fluctuating stresslet. With the inclusion of the fluctuating stresslet, the analytic calculations of the viscosity for an isolated force-free rigid spherical particle under equilibrium and non-equilibrium conditions are now in agreement, both yielding $B=5/2$ (and this result is unaffected by the full derivation below).

The interpretation of the stochastic differential equation of motion, Eq. (7), requires further attention to resolve an ambiguity: the impact of the Brownian force, as determined by the resistance matrix in the fluctuation-dissipation theorem and the mobility matrix in the equation of motion, varies with the coordinates $\mathbf{Q}$ of the particle, while the coordinates change due to this Brownian force. In the Itô interpretation, i.e. using only parameter values at time $t$ just before the impact of the Brownian force, the integration of the equation of motion over a time step $\Delta t$ results in coordinate increments Palanisamy19

[TABLE]

where $g_{Q}$ denotes the metric of the coordinate space, $\bm{\rho}$ is the symmetric tensor solving $\bm{\rho}^{2}=2\pi k_{B}T\bm{\xi}$, and $\bm{\Theta}^{U}(t)$ and $\bm{\Theta}^{E}(t)$ are random vectors with zero mean, unit variance and devoid of correlations, containing six and five unique elements, respectively. The divergence term brings into account the coordinate-dependence of the hydrodynamic matrices Ermak78 ; Kampen ; Ottinger ; Gardiner ; these additional displacements are not evident from Eq. (7) but are crucial to obtaining the equilibrium Boltzmann distribution and, as will be shown below, also contribute to the stress. An alternative interpretation, due to Einstein, gave rise to the name ‘thermodynamic force’ Einstein05 ; Batchelor76 ; GuazzelliMorris .

On a technical note, since in the current context the stress and strain rate tensors are symmetric and traceless $(3\times 3)$ matrices, it proves convenient to replace both by five-vectors so the usual mathematical and numerical techniques can be applied to the resulting symmetric $(11\times 11)$ hydrodynamic matrices KimKarrila ; Palanisamy18 . Because the hydrodynamic matrices and the conservative potential are typically expressed in terms of Cartesian velocities and Cartesian forces, and in angular velocities and torques around Cartesian axes, we furthermore take the freedom of evaluating the r.h.s. of Eqs. (7) and (8) in Cartesian coordinates, henceforth collectively denoted as $\mathbf{X}$. Since the Cartesian angular velocities are not time derivatives of angular coordinates, the rotation angle increments still require transformation to proper generalized coordinates $\mathbf{Q}$ describing the orientation of the colloid in terms of e.g. Euler angles or quaternions, or one may directly update the rotation matrix between the colloid-based axes frame and the space-based axes frame. The latter two options have the advantage that they do not require corrections resulting from the metric. Makino04motion ; Ilie15 ; Palanisamy18 .

Continuing in the Itô representation, we find by some mathematical manipulations of Eqs. (6) through Eq. (8) that the average stress exerted by the colloid on the fluid during a time step $\Delta t$ reads as Palanisamy19

[TABLE]

The first and second term on the r.h.s. are the two deterministic contributions to Eq. (7). The third term, usually derived along another route, results here by combining a term related to the divergence in Eq. (8) with a term resulting from correlations between $\delta\mathbf{F}$ and $\delta\mathbf{S}_{\delta\mathbf{F}}$, and is referred to in the literature as ‘Brownian stress’ or ‘diffusion stress’ Batchelor77 ; Rallison78 ; Bossis89 ; Makino04visco . The fluctuating terms, collected between curly brackets in the last term to Eq. (9), have zero average and may therefore be ignored when calculating the time-averaged stress of a system under shear flow, but their correlations are crucial when applying the Green-Kubo formalism to a quiescent system. In both non-equilibrium and equilibrium cases, the time-averaged correlations of the fluctuating forces and the fluctuating stresslets give rise to two additional non-fluctuating stress contributions, the two $\bm{\nabla}_{\!X}\bm{\rho}^{S}_{\circ}$ terms in Eq. (9), which, to the best of our knowledge, are derived and reported here for the first time.

IV Numerical example

As an illustration of the revised stress calculation, we present numerical simulations of isolated spheroidal particles. For a rigid particle, the hydrodynamic matrices are constant in the body frame and rotate with the body in the space frame, which permits efficient calculations of the motion and the stress without demanding re-evaluation of the hydrodynamics Makino04motion ; Ilie15 ; Palanisamy18 . With $\mathbf{A}$ denoting the body-to-space rotation matrix, the three derivative-containing terms in Eq. (9) turn out to be of the form $\mathbf{A}\mathbf{G}\mathbf{A}^{\mathrm{T}}$, where $\mathbf{G}$ is a body-dependent constant $(3\times 3)$ matrix and the superscript $\mathrm{T}$ indicates transposition. Simulations of a spherical particle confirm the theoretical results mentioned above, with the classical approach of Eq. (3) yielding Einstein coefficients of 5/2 under shear and zero by the Green-Kubo method, while the amendments proposed here yield $B=5/2$ for both non-equilibrium and equilibrium simulations. Numerical results for the Einstein coefficients of spheroidal particles are presented in Fig. 1 as function of the aspect ratio $p=L/D$ between the length $L$ along the symmetry axis and the diameter $D$ in the perpendicular direction, ranging from disk-like oblate to needle-like prolate. For simulations based on Eq. (3), applying the Green-Kubo formalism to quiescent systems yields Einstein coefficients between four ($p\ll 1$) and eight ($p\gg 1$) times higher than the values obtained from the average stress in sheared systems, with an intermediate dip of $B$ approaching zero for near-spherical particles. Simulations based on Eq. (7), however, show consistency between equilibrium and non-equilibrium viscosity evaluations, see Fig. 1. The biggest difference between the classical and the proposed methods is in the Einstein coefficients deduced from the thermal stress fluctuations in quiescent fluids, and mainly results from the inclusion of the Brownian fluctuating stress in Eq. (9). The two novel $\bm{\nabla}_{\!X}\bm{\rho}^{S}_{\circ}$ terms in said expression introduce a relatively modest increase of the Einstein coefficient obtained under shear, by about 1% at $p=20$ and 5% at $p=100$. The impact of these terms is larger for less symmetric bodies, amounting to about $+10\%$ for a semi-disk with diameter-to-thickness ratio of 40 and about $-8\%$ for a helix inscribing 7.5 revolutions in a cylinder with a length-to-diameter ratio of five.

V Conclusions

The perpetual thermal motion of fluids contributes to the viscosity of colloidal suspensions, both by causing the Brownian motion of the colloids and by inducing flucuating stresses on the colloids. Inclusions of these Brownian stresses, absent in current theoretical and numerical implicit solvent methods for suspensions, is therefore necessary to obtain the correct viscosity. Our theoretical analysis of dilute suspensions of spherical particles and numerical simulations of spheroids illustrate the validity of the amendments proposed in the Brownian motion and stress calculation of Eq. (7) and the extended fluctuation-dissipation theorem of Eq. (6). Correlations between the various fluctuating terms then give rise to novel non-fluctuating contributions in the Itô representation of the stress, see Eq. (9). A more detailed exposition of the derivations outlined above, along with additional numerical results on colloids of various shapes, will be presented elsewhere Palanisamy19 . The expression for the time-averaged stresslet on a colloid is readily extended to a collection of $N$ particles, by enlarging the mobility and resistance matrices to $(11N\times 11N)$ matrices including hydrodynamic interactions between all colloids; the total deviatoric stress in the system is then obtained by adding up the stresslets of the individual particles, the virial term due to generalized conservative forces on the colloids, and the shear resistance of the suspending fluid.

Acknowledgments

We thank Prof. Stefan Luding for stimulating discussions. This work is part of the Computational Sciences for Energy Research Industrial Partnership Programme co-financed by Shell Global Solutions B.V. and the Netherlands Organisation for Scientific Research (NWO).