3370
J. Phys. Chem. B 2007, 111, 3370-3378
Cell Models for the Primary Electroviscous Effect J. D. Sherwood† Schlumberger Cambridge Research, High Cross, Madingley Road, Cambridge CB3 0EL, United Kingdom ReceiVed: September 8, 2006; In Final Form: January 12, 2007
The primary electroviscous effect in a nondilute suspension of charged spherical particles is studied by means of cell models. The governing equations are derived, and then analytic results are obtained by restricting attention to the limit of thin double layers, small Hartmann and Peclet numbers, and small potentials. Previous work has assumed that the velocity at the outer boundary of the cell is identical to the imposed flow, as proposed by Simha (J. Appl. Phys. 1952, 23, 1020). Results with this boundary condition are compared against those predicted when the tangential shear stress on the outer boundary is assumed to be unperturbed, as proposed by Happel (J. Appl. Phys. 1957, 28, 1288). Both the hydrodynamic and electroviscous contributions to the effective viscosity are smaller with the Happel boundary condition, showing that such cell models offer a range of predictions and should be used with caution.
1. Introduction Particle-particle interactions play an important role in determining the macroscopic properties of concentrated suspensions of solid particles in liquid. Theories that ignore such interactions are therefore usually valid only in the limit of small solid volume fraction φ , 1. Thus Einstein’s analysis1,2 of the flow field around a single spherical particle in unbounded fluid predicts the viscosity of the suspension correct only to O(φ), and analyses3,4 of the interaction between pairs of particles extend Einstein’s result only to O(φ2). Exact results for higher concentrations are difficult to obtain. Modern computing power sometimes allows direct computation,5 but in the past various effective medium techniques were developed to give some indication of behavior at high concentration.6 One such technique is the cell model, based on a unit cell consisting of a single particle surrounded by a volume of suspending fluid chosen to give the desired solid volume fraction φ. If the particle is a sphere of radius a, the natural choice of cell will be a sphere of radius aφ-1/3, concentric with the particle (Figure 1). Although engineers are glad to get some prediction of the properties of concentrated suspensions, these predictions can be strongly affected by the choice of boundary conditions on the outer surface of the cell, and must be treated with caution. It is therefore useful to investigate a range of boundary conditions in order to have some feeling for the uncertainty of the predictions. Our interest here lies in the viscosity of a concentrated suspension of spheres, and in the contribution of ionic double layers to this viscosity. We restrict our attention to the primary electroviscous effect, caused by deformation of the charge cloud around a single spherical particle.7 If the electrical double layers of two particles overlap, there will be a (repulsive) force between the particles even in the absence of flow, and this force strongly modifies the distribution of particles and the interparticle stresses (the secondary electroviscous effect). To avoid such effects, we assume that the Debye length κ-1, which characterizes the dimension of the charge cloud, is small compared to the interparticle separation O(aφ-1/3). †
E-mail:
[email protected].
Figure 1. The cell of radius aR concentric with the spherical particle of radius a.
Cell models for the viscosity of suspension of uncharged spheres were investigated by Simha8 and by Happel:9 their results are discussed by Sherwood10 and reviewed in section 4.3. The fluid velocity around the sphere at the center of the cell is modified from that around a sphere in unbounded fluid, and in consequence the additional stresses caused by the deformation of an ionic charge cloud around a charged particle at the center of the cell differ from the stresses in unbounded fluid. This was investigated by Ruiz-Reina et al.11 The electrical stress depends upon the dimensions of the cell, and we interpret this in terms of an electrical contribution to the suspension viscosity that depends (nonlinearly) upon the solid volume fraction φ. In a dilute suspension the primary electroviscous effect gives rise to an O(φ) change in the mean stress, and the secondary electroviscous effect (electrical interactions between particles) gives an O(φ2) contribution. We shall nevertheless follow the original classification of Dobry7 and classify stresses due to the hydrodynamic deformation of the charge cloud around a single particle as part of the primary electroviscous effect, whether linear in φ or not. Our aim here is to show how these stresses are modified by the choice of boundary condition at the outer spherical boundary of the cell. It should be noted that the use of cell models to predict electrophoresis (rather than the electroviscous effect) in con-
10.1021/jp065862m CCC: $37.00 © 2007 American Chemical Society Published on Web 03/14/2007
Cell Models for the Primary Electroviscous Effect
J. Phys. Chem. B, Vol. 111, No. 13, 2007 3371
centrated suspensions has a long history:12-19 recent literature is reviewed by Carrique et al.18 and the choice of boundary condition on the outer boundary of the cell is discussed by Ahualli et al.19 Most authors adopt either Happel’s condition of zero tangential shear stress on the outer boundary20 or Kuwabara’s zero vorticity condition.21 The difference between predictions is small when the double layer is thin, in which case vorticity outside the double layer is zero and Kuwabara’s zero vorticity boundary condition gives the exact result for an isolated particle.12 More generally, when the Debye length κ-1 is comparable with the particle radius a, vorticity is no longer zero away from the particle, and it is less obvious that Kuwabara’s boundary condition is correct.17 Various authors13,17 have discussed the choice between Dirichelet or Neumann boundary conditions for the electrochemical potential of ions at the outer boundary of the cell in the context of electrophoresis. In the electroviscous problem considered here we adopt a condition of zero derivative normal to the outer boundary for the electrochemical potential. Electrochemically driven fluxes of ions into or out of the cell are zero: this is intended to represent the absence of large-scale gradients of electrical potential or ionic concentration within the suspension. In section 2 we summarize the equations governing the equilibrium charge cloud, and in section 3 we consider the deformation of the charge cloud caused by fluid flow. Section 4 determines the average stress in the suspension. In the limit aκ . 1 it is possible to obtain asymptotic results analytically using methods first developed by Dukhin and Shilov22 and Fixman:23 the primary electroviscous effect around a particle in unbounded fluid was studied in this limit by Hinch and Sherwood.24 In section 5 this analysis is applied to a charged particle within a cell, and the electroviscous effect is determined for the two different sets of boundary conditions on the outer boundary of the cell.
the Poisson-Boltzmann equation becomes I
∇ p)2
∑ i)1
κ-1 )
( ) kT
∑i e
∇ φ ) -F/
∂φ0 )0 | ∂r r)aR
ni0
)
ni∞
(
)
eziφ0 exp kT
(3)
is the number density of the ith species far from any where charged surfaces, φ0 is the equilibrium potential, and kT is the Boltzmann temperature. Inserting the Boltzmann distribution (3) into the Poisson equation (2) leads to the PoissonBoltzmann equation for the equilibrium potential φ0. In terms of the nondimensional equilibrium potential
p)
eφ0 kT
A + exp(-η) A - exp(-η)
(4)
(9)
where η ) (r - a)κ is a nondimensional distance from the charged surface, which we take to be at the slipping plane where φ ) ζ. This boundary condition implies that
A)
exp(h) + 1 exp(h) - 1
(10)
eζ 2kT
(11)
where
3. The Deformed Charge Cloud 3.1. General Equations. The ions move under the influence of electric and thermodynamic forces with velocity
vi ) u + ωi(-ezi∇φ - kT∇ ln ni) ) u - ωi∇ψI
ni∞
(8)
However, as discussed above, we require that there should be no overlap of charge clouds, i.e., we assume κ-1 , a(R - 1), and there will be little difference between (8) and the simpler condition φ0 f 0 as r f ∞. We shall later assume that the Debye length κ-1 is small compared to the particle radius a. In this limit, to a first approximation we may neglect curvature and use analytic solutions of the plane one-dimensional Poisson-Boltzmann for symmetric electrolytes:
(2)
where is the permittivity of the electrolyte. The equilibrium ionic number density ni0 is given by the Boltzmann distribution
(7)
The potential φ0 around an isolated particle in unbounded fluid decays as r-1 exp(-κr) as the distance r from the center of the particle becomes large. In a bounded cell of radius aR ) aφ-1/3 the natural boundary condition to impose is
h)
2
(6)
(zi)2ni∞
∇2p ) -κ2 sinh p
(1)
where -e is the electron charge. The electric potential φ satisfies the Poisson equation
1/2
For a monovalent symmetrical electrolyte, with z(1) ) -z(2) ) 1, the Poisson-Boltzmann equation (5) becomes
I
nizie ∑ i)1
(5)
2
exp(p/2) )
We assume that the suspending electrolyte contains I species of ions, each with number density ni and valence zi. The charge density F is therefore
exp(-zip)
kT
The equilibrium charge cloud thickness is characterized by the Debye length
2. The Equilibrium Charge Cloud
F)
( ) e2zini∞
(12a) (12b)
where u is the fluid velocity, ωi is the mobility of the ith species of ion within the (dilute) electrolyte, and
ψi ) eziφ + kT ln ni
(13)
is the electrochemical potential of the ith species of ion (at constant pressure). We assume that the ions are not taking part in any reactions, so that in steady state the ionic number densities satisfy the conservation equation
3372 J. Phys. Chem. B, Vol. 111, No. 13, 2007
Sherwood
∇‚(niu - ωi(ezini∇φ + kT∇ni)) ) 0
∇‚u ) 0
(14)
(23)
If U is a typical fluid velocity and ω a typical ion mobility, the Pe´clet number
Colloidal particles are sufficiently small that inertia is negligible, and so u satisfies the Stokes equation
Pe ) U/(ωkTκ)
0 ) -∇P + µ0∇2u - F∇φ
(15)
measures the ratio of the convective forces on the ions to the diffusive forces. If Pe , 1 (as is usually the case) the charge cloud is only slightly deformed from equilibrium. We can therefore expand the number densities and potential as
(24)
where P is the fluid pressure and the last term on the righthand side of (24) represents the electrical forces acting on the fluid. Expanding F∇φ in powers of the Peclet number, the Stokes eq 24 becomes
ni ) ni0 + ni1 + ...
(16a)
φ ) φ0 + φ1 + ...
(16b)
2) - en∞φ1(ep - e-p)] (25)
where ni1 and φ1 are O(Pe). We restrict ourselves to a 1-1 electrolyte, and use superscripts ( to denote cations and anions. The average mobility is ω ) (1/2)(ω+ + ω-), though later we shall set ω+ ) ω- ) ω to simplify the analysis. The gradient of the ionic electrochemical potentials ψ( may now be expanded as
and various terms on the right-hand side can be absorbed into the pressure P. 3.2. Equations for a Straining Flow. We now restrict our attention to a spherical particle, and nondimensionalize all lengths by the particle radius a, so that the outer boundary of the cell is at r ) R. The nondimensional volume of the cell is
∇ψ( ) (e∇(φ0 + φ1) +
)
( kT∇(n( 0 + n1 )
n( 0
n( 1
+
4 V ) πR3 3
+ O(Pe)2
kT ( ∇χ + O(Pe)2 n∞
(17a)
( )
χ ( ) n( 1 exp((p) ( φ1
en∞ kT
φ ) R-3
ω(kT ( ∇χ n∞
(19)
To leading order in the Pe´clet number the ionic fluxes are
n(v( ) exp(-p)(n∞u - ω(kT∇χ()
(20)
∇2χ- + ∇p‚∇χ+ ) +
n∞ ω+kT n∞ ω-kT
u‚∇p
(21a)
u‚∇p
(21b)
There can be fluxes of ions into or out of any Stern layer on the particle surface, and the effect of these fluxes on the electroviscous effect has been examined both for a particle in unbounded fluid26,27 and for a cell model.28 However, we assume here that no such layer is present. If n is the normal to the particle surface, the normal velocity n‚u ) 0 at the surface and hence zero flux of ions into the surface of the particle requires
n‚v( ) n‚∇χ( ) 0
u ) E‚x
(28)
in the absence of any solid particles. This is modified by the presence of the particle to produce a velocity field of the form
x‚E‚x f(r) r2
u ) E‚xg(r) + x
(29)
Incompressibility (23) implies
g′ + f′ + 3fr-1 ) 0
(30)
where we use the notation ′ ≡ d/dr. The corresponding form for the electrochemical potentials is a quadrupole
and the ion conservation (14) becomes
∇2χ+ - ∇p‚∇χ- ) -
(27)
Any vorticity in the imposed flow field merely rotates the particle and surrounding fluid, so we consider a straining flow that would be
(18)
The velocities v( of the ions (12) can therefore be written as
(26)
and the solid volume fraction in the cell is
(17b)
where the zeroth order (equilibrium) terms in (17a) have vanished, by (5), and as in ref 25 we use the notation
v( ) u -
µ0∇2u ) kT(χ+e-p - χ-ep)∇p + ∇[P - n∞kT(ep + e-p -
(22)
The fluid is assumed incompressible, so that the fluid velocity u satisfies the equation of continuity
χ( )
( )
n∞a2 x‚E‚x q((r) ω(kT r2
(31)
where, by (21), the q( satisfy
q′′( + 2r-1q′′( - 6r-2q( - p′q′′( ) - rp′(g + f)
(32)
Outside the equilibrium charge cloud, where p ≈ p′ ≈ 0, the right-hand side of (32) is small, and solutions of the homogeneous equation have the form
q( )
c5 r3
+ c6r2
(33)
where c5, c6 are unknown constants that must be determined as part of the solution. The equation of motion (eq 25) becomes
Cell Models for the Primary Electroviscous Effect
f ′′′ +
(
J. Phys. Chem. B, Vol. 111, No. 13, 2007 3373
)
-p 2 9f ′′ 6f ′ 30f kT(aκ) ωq+e ωq_ep dp + 2 - 3 ) 2 2 r r r r µ0e ω ω+ ω- dr (34)
and we note that on the right-hand side of (34) the Hartmann number
H ) kT/(µ0e2ω)
f(r) ) c1r2 +
r3
+
c3 r5
5c1r2 2c3 g(r) ) - 5 + c4 2 5r
(36a)
(36b)
where the ci are unknown constants that must be determined. The pressure corresponding to (36) is
x‚E‚x r2
p ) µ0h(r) where
(
(37)
)
21c1r2 h(r) ) 3 2 r 2c2
(38)
(39)
The boundary condition (22) that there should be zero ionic flux into the surface of the particle requires
q′+(1) ) q′_(1) ) 0
g(R) ) -
c2 R
+ 3
c3 )0 R5
5c1R2 2c3 + c4 ) 1 2 5R5
(41a)
(41b)
Happel9 took one boundary condition at r ) R to be zero perturbation of the fluid velocity normal to the outer boundary, i.e., x‚u - x‚E‚x ) 0, and hence
g(R) + f(R) ) -
3c1R2 c2 3c3 + 3 + 5 + c4 ) 1 2 R 5R
(
)
x‚E‚x x‚e‚x ) E‚x - x 2 2 R R
e‚x‚ - x
r)R
(44)
and this supplies a second boundary condition at r ) R:
c2 8c3 Rg′(R) + g(R) + f(R) ) -4c1R2 + 3 + 5 + c4 ) 1 2 R 5R (45) The governing equations involve only the perturbed electrochemical potentials χ(, rather than the separate perturbations to the ionic number densities n1 and electrical potential φ1. We therefore seek a boundary condition for the O(Pe) electrochemical potentials χ (, which are quadrupoles that decay more slowly than the equilibrium potential φ0 when r . 1. One option might be to assume zero perturbation χ( ) 0 on the outer boundary r ) R (cf. Simha8). However, given the zero gradient condition (8) on the equilibrium potential φ0, it seems more natural to assume n‚∇χ( ) 0, and hence
(46)
This implies, by (19), that ions can still be convected into or out of the cell by the fluid, but electrochemically driven fluxes are zero, as one might expect in a suspension with no mean gradient of electical potential or ionic concentration. Ruiz-Reina and co-workers11,28,29 adopted a mixed boundary condition, with n( 1 ) n‚∇φ1 ) 0 on r ) R. The physical significance of this mixed condition is not clear, but physical intuition suggests that it is more constraining than (46), and so leads to a larger estimate of the electroviscous effect. 4. The Suspension Viscosity
(40)
However, the appropriate boundary conditions on the outer surface r ) R of the cell are less evident. Simha8 assumed that the fluid velocity on the outer surface of the cell is unperturbed, and given by (28), and hence
f(R) ) c1R2 +
)
g′ 2f + [E‚xx + xE‚x] + 2gE + r r2 4f x‚E‚x 2f 2f ′ - xx 3 + 2Ix‚E‚x (43) r r r
q′(R) ) 0
3.3. Boundary Conditions. Equations 21 and 25 for χ( and u require boundary conditions. The particle is at rest, and to ensure u ) 0 at its surface r ) 1 we require
f(1) ) g(1) ) 0
(
The perturbation to the tangential shear stress over the outer boundary S1 of the cell is zero if
(35)
characterizes the ratio of electrical forces acting on the charge cloud to viscous forces. Outside the equilibrium charge cloud, solutions of (30) and (34) for the fluid velocity have the form
c2
2e ) ∇u + (∇u)T )
(42)
For his second boundary condition Happel assumed that there is zero perturbation to the tangential shear stress over the outer surface of the cell. The rate of strain e of the perturbed velocity (29) is
4.1. The Average Strain Rate. The average strain rate within the cell is given by the volume average
2〈eij〉 )
1 V
(
∂u
∂u
)
∫V ∂xji + ∂xij
dV )
∫
1 (u dS + uj dSi) (47) V S i j
where V ) (4/3)πR3 is the volume occupied by the entire cell, with outer boundary S. Taking dSi ) R2xˆ i dΩ, where dΩ is a unit of solid angle, and using the results
δ, ∫xˆ ixˆ j dΩ ) 4π 3 ij (δ δ + δikδjl + δilδjk) ∫xˆ ixˆ jxˆ kxˆ l dΩ ) 4π 15 ij kl
(48)
we find from (47)
2〈eij〉 )
3 2π
∫V(g(R)Einxˆ nxˆ j +
(
f(R)xˆ ixˆ jxˆ mEmnxˆ n) dΩ ) 2Eij g(R) +
(
) 2Eij c4 -
)
21c1R2 2c2 + 3 10 5R
)
2f(R) 5
(49)
3374 J. Phys. Chem. B, Vol. 111, No. 13, 2007
Sherwood
4.2. The Average Deviatoric Stress. The total stress is σ + m, where the Maxwell stress tensor
1 mij ) EiEj - E2δij 2
(50)
and in the absence of inertia ∇‚(σ + m) ) 0. Hence the volume averaged deviatoric stress is
〈σij + mij〉 )
Ω1(R) ) 4R10 - 25R7 + 42R5 - 25R3 + 4
(51)
{(
3µ0 4π
∫S
ΦS )
)
}
µ0Eij {7g′(R)R + 10f(R) + 10g(R) + 4f ′(R)R - 2h(R)} 5
(
6c2 µ0Eij -21c1R2 - 3 + 10c4 ) 5 R
)
c1 2R3 , c2 ) - (5R7 + 2) Ω2 2
(59a)
c3 )
5c1R7 c1 , c4 ) (2R7 + 5) 2 2
(59b)
Ω2(R) ) 2R10 - 5R7 + 5R3 - 2
(60)
where
and we obtain the effective viscosity
µ ) 1 + ΦH µ0 where
]
[
]
10c4R - 21c1R + 4c2 5
f(1) ) c1 + c2 + c3 ) 0
(53)
(62)
(
5c1 2c3 + c4 ) 0 2 5
(54a) (54b)
10R3(R2 - 1) 1 - R7 , c2 ) c 1 2 Ω1 R -1
(55a)
(63)
Both (57) and (61) tend to Einstein’s viscosity
5 µ ) µ0 1 + φ + O(φ2) 2
If we combine (54) with Simha’s boundary conditions (41) on the outer boundary of the cell, we find
c1 )
(25R7 + 10) Ω3
Ω3(R) ) 10R10 - 10R7 - 21R5 + 25R3 - 4
4.3. Results for Uncharged Particles. In the absence of any charge cloud, the forms (36) for f and g hold everywhere within the fluid, and so the boundary condition u ) 0 on the surface of the particle implies
g(1) ) -
ΦH ) with
10c2 3
(61)
(52)
-21c1R5 - 6c2 + 10c4R3 〈σ〉 ) ) µ0 µ) 2〈e〉 10c4R3 - 21c1R5 + 4c2 µ0 1 -
(58)
c1 )
Both the volume-averaged deviatoric stress (52) and strain rate (49) are linear in Eij, so that we can define the effective viscosity µ of the suspension as
[
(10R7 - 10) Ω1
This is Simha’s8 result. If we instead adopt Happel’s boundary conditions (42) and (45) on the outer boundary of the cell, we find
)
(
(57)
where
g′ 2f [Eimxmxk + + r r2
xmEmnxn 4f xixk + xiEkmxm] + 2gEik + 2f ′ r r3 xmEmnxn (2f - h)δik xˆ jxˆ k dΩ r2
(56)
and the effective viscosity is
µ ) 1 + ΦS µ0
However, the expansion (16b) implies that the O(Pe) contribution to the Maxwell stress consists of terms of the form ∇φ0∇φ1 that have decayed exponentially at r ) R, and can be ignored since we have assumed (aκ)-1 , R - 1. An analysis that does not make this assumption has been presented by Ruiz-Reina et al.29 The hydrodynamic stress in (51) is σ ) 2µ0e - PI, with the strain-rate e outside the charge cloud given by (43) and the pressure P given by (37). Hence
)
where
3 4πR3
3 (σ + mik)xj dSk ∫V∂x∂ k[(σik + mik)xj] dV ) 4πR 3∫S ik
〈σij + mij〉 )
()
c1 4R7 + 21R2 - 25 R7 - R 2 , c4 ) (55b) 2 10 R -1 R2 - 1
c3 ) c1
)
φ,1
(64)
in the low volume fraction limit φ , 1, correcting Happel’s9 result. Both models predict a viscosity that becomes infinite in the limit R f 1, even though a suspension of monodisperse spheres will have jammed solid well before this limit is attained. Happel’s constraint on the shear stress at r ) R is weaker than Simha’s constraint on the fluid velocity, and consequently Happel’s viscosity (eq 61) is smaller than that of Simha (eq 57). The predicted increases in viscosity ΦS (eq 57) and ΦH (eq 61) are shown in Figure 2, and are discussed further by Sherwood.10 Note that we have chosen to follow Batchelor30 and compute the volume averages of stress 〈σij〉 (51) and strain rate 〈eij〉 (47). The viscosity (53) is then defined as the ratio 〈σij〉/(2〈eij〉). An alternative route1,9,31 is to compute the rate of energy dissipation Φ. This is given by Batchelor30 (his eq 4.8). In the absence of
Cell Models for the Primary Electroviscous Effect
J. Phys. Chem. B, Vol. 111, No. 13, 2007 3375
Figure 2. Contributions to the effective viscosity, as a function of solid volume fraction φ. Hydrodynamic contribution: (a) ΦH (Happel boundary condition) and (b) ΦS (Simha boundary condition). Electroviscous contribution: (c) ΦHe (Happel) and (d) ΦSe (Simha).
electrical effects, inertia, body forces, and couples, the stress is symmetric, and
1 Φ) V
∂ui ∂〈ui〉 1 〈σij〉 + σ dV ) V ∂x ij ∂xj V j
∫
) 〈eij〉〈σij〉 +
∂(u′iσij) dV (65) V ∂xj
∫
∂〈u′iσij〉 ∂xj
(66)
where u′ is a fluctuation with zero mean. If the average is taken over a volume V such that the suspension is statistically uniform, the term ∂j〈u′iσij〉 in (66) vanishes,30 and we conclude that the mean rate of energy dissipation is given by the product of the mean stress and mean strain rate. However, the sample of the suspension inside the spherical cell is not statistically uniform. Integrating (65) by means of the divergence theorem gives
Φ ) 〈eij〉〈σij〉 +
∫
1 u′σ dS V S i ij j
(67)
If the velocity fluctuation
(
x‚E‚x 2 - E‚x 5 r2
u′ ) u - 〈e〉‚x ) f(r) x
)
(68)
6µ0 1 u′‚σ dS ) f(R)(Rg′(R) + 2Rf ′(R) - h(R)) (69) V S 25
∫
The energy analysis9,31 neglects the surface integral on the righthand side of (67), and assumes that Φ ) 〈e〉:〈σ〉 ) 2µ〈e〉:〈e〉, i.e., it assumes that the average of the product of stress and strain rate is equal to the product of their averages. The average viscosity is then taken to be
(70)
In my earlier work10 on suspensions of uncharged spheres I noted that Happel had incorrectly assumed 〈eij〉 ) Eij, whereas (49) with the Happel boundary conditions gives
3 〈eij〉 ≈ Eij 1 + φ + ... for φ , 1 2
(
)
5. Thin Double Layer (aK)-1 , R - 1 5.1. The Velocity Field Outside the Equilibrium Charge Cloud. We follow an analysis very similar to that of Hinch and Sherwood.24 Outside the double layer the perturbed quantities take the form of quadrupoles / n( 1 ) ni
(
(71)
This modified strain rate (eq 71) did not suffice to explain the discrepancy between Happel’s viscosity and that of Einstein (eq 64) at low volume fractions when the viscosity was evaluated as the ratio (eq 53), and I erroneously suggested the possibility of an algebraic error in Happel’s analysis. However, I neglected that Happel used an energy formulation, an error pointed out
)
(
)
3r2 3r2 x‚E‚x 1 x‚E‚.x 1 + 5 , φ1 ) φ/1 2 + 5 2 3 3 r r 2R r r 2R
(72)
so that the boundary condition (46) is satisfied at r ) R, with ∑izin/i ) 0 to ensure zero charge density outside the double layer. Derivatives normal to the surface of the particle across the thin charge cloud of thickness κ-1 are typically O(aκ) larger than tangential derivatives, so that tangential fluxes of ions are sufficiently small that they do not affect the radial flux balance at leading order. The perturbed ion densities and potential are therefore as in ref 24 with
(
ni1 ) exp -
is zero over the surface of the cell, as in Simha’s analysis,8 the surface integral is zero, but in general
µ ) Φ/(2〈e〉:〈e〉)
to me by an anonymous referee. Fluctuations in strain rate and stress in the integral on the right-hand side of (67) are both O(φ), so that this integral is O(φ2). Consequently Happel’s Φ ) 2µEijEji(1 + 5.5φ + ...), combined with (70) and (71), leads to Einstein’s result (eq 64) at O(φ), as found by Zholkovskiy et al.31 The integral on the right-hand side of (67) cannot be neglected at higher solid volume fractions. This leads to problems with the energy formulation, discussed by Zholkovskiy et al.31 We could now solve eqs 32 and 34 numerically, for a range of ζ, aκ, H, and R, and investigate the effect of the choice of boundary condition on the outer surface r ) R of the cell. This would be a considerable undertaking if done exhaustively. However, much insight can be gained by studying the important (and analytically tractable) limit of a thin double layer (aκ f ∞), the approach adopted in the next section.
)(
)
eziφ0 i ezi n1(+) - ni∞(φ1 - φ1(+)) kT kT
(73)
i where n1(+) and φ1(+) are values just outside the double layer:
(
i n1(+) ) n/i x‚E‚x 1 +
)
(
)
3 3 , φ1(+) ) φ/1x‚E‚x 1 + 5 2R5 2R (74)
The perturbed ion density (73) represents a Boltzmann distribution with the perturbed potential φ0 + φ1. The electrochemical potential does not vary across the double layer (to leading order). We restrict our attention to a 1-1 electrolyte, so that n/2 ) n/1. The fluid velocity differs from that around an uncharged sphere because of the electrical forces acting within the perturbed charge cloud. The velocity immediately adjacent to the charged surface can be obtained by integrating the Stokes eq 25 within the (thin) charge cloud. Details are given by Hinch and Sherwood,24 who compute the electroosmotic contribution to the tangential velocity utan in terms of a slip velocity of the form
(
)
x‚E‚x Vs(r) r2
E‚x - x
(75)
The electroosmotic tangential velocity Vs tends to an apparent slip velocity V∞s outside the charge cloud. However, the presence of an outer boundary at r ) R has modified the
3376 J. Phys. Chem. B, Vol. 111, No. 13, 2007
Sherwood
electrochemical potential outside the charge cloud by a factor of 1 + (3/2)R-5 (eq 74). Consequently, the tangential velocity given by Hinch and Sherwood24 is modified by this same factor, and the apparent slip velocity immediately outside the charge cloud becomes
V∞s
[
8 A + (kTn/1 + )(kTn/1 - en∞φ/1) ln 2 A 1 µ0κ A 3 en∞φ/1) ln 1 + 5 (76) A+1 2R
)(
)
Outside the double layer the fluid velocity is given by (29), with f and g as in (36). To leading order (since the charge cloud is thin) the normal velocity u‚x ) 0, so that
g(1) + f(1) ) -
3c3 3c1 + c2 + + c4 ) 0 2 5
(
(78)
The boundary condition (78) implies that c3 and c1 are modified by the nonzero slip velocity Vs. However, r - 1 is O(aκ)-1 within the charge cloud, and so in (79) O(Vs) contributions to 2c3 - 5c1 can be neglected compared to the term Vs. The charge cloud does not affect the boundary conditions at the outer boundary r ) R of the cell. If we take the Simha outer boundary conditions (41) at r ) R, together with (77) and (78) we obtain
(80a) (80b)
Ω1c3 ) 10(R10 - R5) - V∞s (10R10 - 25R7 - 15R5) (80c) Ω1c4 ) 4R10 + 21R5 - 25R3 + V∞s (-15R7 + 21R5 - 6) (80d) where Ω1(R) is defined by (56). The effective viscosity (53) becomes 7 2 ∞ 10(R7 - 1) Vs (6R - 21R + 15) µ )1+ µ0 Ω1 Ω1
where Ω2(R) is defined by (60). The effective viscosity (53) becomes µ )1+ µ0
25R7 + 10 - 15V∞s (R7 - 1)
10R10 - 10R7 - 21R5 + 25R3 - 4 - V∞s (9R7 - 21R5 + 21R2 - 9)
≈1+
∞ 25R7 + 10 15Vs Ω4 Ω3 (Ω )2
(83)
Ω4 ) 10R17 - 25R14 + 14R12 + 15R10 - 35R9 + 15R7 + 35R5 - 25R3 - 14R2 + 10 (84)
)
Ω1c2 ) 10(R3 - R10) + V∞s (6R10 - 21R5 + 15R3)
(82d)
where Ω3(R) is defined by (63) and
xu‚x x‚E‚x ) E‚x - x 2 [Vs + (2c3 - 5c1)(r 2 r r 1) + ...] r - 1 , 1 (79)
Ω1c1 ) 10(R5 - R3) + V∞s (10R3 - 4 - 6R5)
Ω2c4 ) 2R10 + 5R3 + V∞s (3 - 3R7)
(77)
Combining the electroosmotically generated tangential velocity Vs(r) with the inner expansion of the velocity (29) outside the double layer, we find that the tangential velocity within the double layer takes the form
utan ) u -
(82c)
3
However, for the tangential fluid velocity immediately outside the double layer to match onto the slip velocity V∞s we require
5c1 2c3 V∞s ) g(1) ) + c4 2 5
Ω2c3 ) 5R10 - 5V∞s (R10 - R7)
and where we have assumed V∞s , 1. We now know the effective viscosity of the cell (53) in terms of the slip velocity V∞s , and by (76) we know V∞s in terms of the coefficients n/1, φ/1 of the quadrupoles (72) outside the double layer. It remains to determine these coefficients. 5.2. The Flux Balance Within the Deformed Double Layer. The constants n/i and φ* are determined by considering fluxes of ions within the double layer. Integrating the ion conservation eq 14 across the thickness of the double layer, we obtain the flux balance
∫dl∇‚(uni0) dr ) ∫dl∇t‚ωi(ezini0∇tφ1 +
(
kT∇tni1) dr + ωi ezini∞
|
∂φ1 ∂r
| )
∂ni1 ∂r
dl+ + kT
dl+
(85)
where dl stands for the double layer, and dl+ for just outside the double layer. On the right-hand side of (85) the first term represents the divergence of the tangential flow of ions within the charge cloud due to the tangential gradient of electrochemical potential. An explicit expression for this is evaluated by Hinch and Sherwood,24 but their result needs to be increased by a factor 1 + (3/2)R-5 because of the modified electrochemical potential (74) just outside the double layer. The second term on the right-hand side of (85) represents the flux of ions into the double layer due to gradients of electrochemical potential normal to the surface of the particle, and is
xˆ ‚∇(kTn( 1 ( en∞φ1) ) -3(1 - R-5)ω((kTn/1 ( en∞φ/1)xˆ ‚E‚x
r ) 1 (86)
(81)
Alternatively, for the Happel boundary conditions (42) and (45) at r ) R
Ω2c1 ) 2R3 - 2V∞s (R3 - 1)
(82a)
Ω2c2 ) -5R10 - 2R3 + V∞s (3R10 - 3R3)
(82b)
which differs from the corresponding result of Hinch and Sherwood24 by the factor 1 - R-5. The left-hand side of (85) represents the divergence of tangential convection of ions, with the tangential velocity given by (79). The electroosmotically generated slip velocity (76) differs from that of Hinch and Sherwood24 by a factor 1 + (3/ 2)R-5, whereas the term in r - 1 in (79), required by matching with the fluid velocity outside the double layer, differs due
Cell Models for the Primary Electroviscous Effect
J. Phys. Chem. B, Vol. 111, No. 13, 2007 3377
to the modified coefficients ci for the external flow given by (80) for the Simha boundary condition or (82) for the Happel condition. We now assume that the particle is positively charged and immersed in a 1-1 electrolyte. Using the results of ref 24 the flux balance (85) for the negative counterions is
-
(
[
)
4ωkTn/1H h 3 A A 1+ 5 e - 1 - 3 ln - ln n∞ A-1 A+1 2R
3 1+ 5 2R
)]
(
(
(aκ) µ0κ ω 1 - R
[
(
)
)
12ω3 ) - 2 (kTn/1 - en∞φ/1)(eh - 1) 1 + 5 aκ 2R 3ω(kTn/1 - en∞φ/1)(1 - R-5) (87) a
)
1+
-5
)(
1+
)
3 2R5
(92)
40R10 - 140R5 + 100R3 + O(Vs) (93) Ω1
8h2n∞ µ ΦSe ) 1 + ΦS + µ0 (aκ)2µ0κ2ω
(94)
where ΦS is given by (58) and
ΦSe )
(40R10 - 140R5 + 100R3)(6R7 - 21R2 + 15) (1 - R-5)Ω12
(
1+
)
3 (95) 2R5
In the dilute limit φ ) R-3 , 1 the viscosity (94) can be expanded as
/
4ωkTn1H h A 3 A - ln e - 1 - 3 ln 5 n∞ A+1 A-1 2R
(
2
and the effective viscosity (81) is
6n∞ (4c3 - 10c1) A ln - 4ωeφ/1H(h + e-h - 1) 2 aκ A + 1 aκ 1+
2
For Simha’s boundary condition, by (80)
where the Hartmann number H ) kT/µ0e2ω. The flux balance for the positive co-ions is
-
(
4c3 - 10c1
8h2n∞
4c3 - 10c1 )
6n∞ (4c3 - 10c1) A ln - 4ωeφ/1H(h - eh + 1) 2 aκ A 1 aκ
(
V∞s ) -
)]
3 2R5
(
)
+
)-
3ω (kTn/1 + en∞φ/1)(1 - R-5) (88) a
If the Hartmann number H is small we may neglect perturbations to the liquid velocity when computing the contribution of convection in the ion conservation equation. A typical value might be H ) 0.29;27 the errors caused by neglecting terms in H in the flux balance have been quantified when the sphere is in unbounded fluid.27 Surface potentials are typically of order |ζ| < 100 mV (i.e., h < 2); here we assume that the equilibrium potential is small, i.e., h , 1. The errors associated with this assumption have been investigated by Sherwood25 and Watterson and White.32 Once we have restricted ouselves to thin double layers and low potentials, the electroviscous effect will inevitably be small. However, this simplification makes it much easier to demonstrate the effect of changing the boundary conditions on the outer surface of the spherical cell. The above flux balances then simplify to become, for the negative counterions
2n∞(4c3 - 10c1) a2κ2
A ) (kTn/1 - en∞φ/1)ω-(1 - R-5) ln A-1 (89)
and for the co-ions
2n∞(4c3 - 10c1) 2 2
aκ
A ln ) (kTn/1 + en∞φ/1)ω+(1 - R-5) A+1 (90)
where, by (10),
1 h A A ) ln (1 + exp(h)) ≈ , ln ) ln A-1 2 2 A+1 h 1 ln (1 + exp(-h)) ≈ - (91) 2 2 For typical 1-1 salts we introduce only a little error if we set ω+ ) ω- ) ω within the electrolyte. We find, by (76)
8h2n∞ µ 125 2 5 375 2 φ + φ + ... )1+ φ+ 15φ + 2 2 µ0 2 8 2 (aκ) µ0κ ω (96)
(
)
which agrees with Booth’s O(φ) result33 for a dilute suspension. For Happel’s boundary condition, by (82)
4c3 - 10c1 )
20R10 - 20R3 + O(Vs) Ω2
(97)
and the effective viscosity (83) is
8h2n∞ µ ΦHe ≈ 1 + ΦH + µ0 (aκ)2µ0κ2ω
(98)
where ΦH is given by (62) and
ΦHe )
(
300(R10 - R3)Ω4 (1 - R
-5
)Ω2Ω23
1+
3 2R5
)
(99)
In the dilute limit φ ) R-3 , 1 the viscosity (98) can be expanded as
8h2n∞ µ 5 5 135 2 φ + ... ) 1 + φ + φ2 + 15φ + µ0 2 2 2 (aκ)2µ0κ2ω (100)
(
)
and we again recover Booth’s result33 for a dilute suspension. In the absence of any electroviscous effect, Happel predicts a lower viscosity than Simha. Figure 2 shows ΦH, ΦS, ΦHe, and ΦSe as functions of the solid volume fraction φ ) R-3. Happel’s viscosity (61) of uncharged spheres is less than that of Simha (57), i.e. ΦH < ΦS, and similarly the electroviscous contribution ΦHe predicted with a Happel boundary condition is less than the electroviscous contribution ΦSe predicted with Simha’s boundary condition. Figure 3 shows the ratios (a) ΦHe/ΦH and (b) ΦSe/ΦS. The electroviscous contributions grow more rapidly than the hydrodynamic contributions, though the electroviscous contributions
3378 J. Phys. Chem. B, Vol. 111, No. 13, 2007
Figure 3. Ratio of electroviscous to hydrodynamic contributions to the effective viscosity, as a function of solid volume fraction φ: (a) ΦHe/ΦH (Happel boundary condition) and (b) ΦSe/ΦS (Simha boundary condition).
will always be small in this limit of thin double layers. We note that the ratio of electroviscous to hydrodynamic contributions is not strongly affected by the choice of boundary condition. 6. Concluding Remarks As far as I know, there are no exact results against which we can compare the electroviscous contributions to the viscosity predicted by cell models. However, various O(φ2) results are available for suspensions of uncharged spheres. In extensional flow Batchelor and Green3 found µ ) µ0(1 + 2.5φ + kφ2 + O(φ3)), with k ) 7.6 in the absence of Brownian motion. At zero Peclet number (when Brownian motion dominates) Batchelor4 found k ) 6.2 and Bergenholtz et al.5 found k ) 5.92. The self-consistent estimate6 µ ) µ0(1 - (5/2)φ)-1 corresponds to k ) 6.25 and experiments34 suggest k ) 4 ( 2. These results lie between the predictions of the cell model with the two boundary conditions. Given that the ratio of electroviscous to hydrodynamic contributions (Figure 3) does not vary markedly with the choice of boundary condition, it seems reasonable to suggest that the correct O(φ2) electroviscous contribution lies somewhere between the predictions using the two different boundary conditions. References and Notes (1) Einstein, A. Ann. Phys. (Weinheim, Ger.) 1906, 19, 289. (2) Einstein, A. Ann. Phys. (Weinheim, Ger.) 1911, 34, 591.
Sherwood (3) Batchelor, G. K.; Green, J. T. J. Fluid Mech. 1972, 56, 401. (4) Batchelor, G. K. J. Fluid Mech. 1977, 83, 97. (5) Bergenholtz, J.; Brady, J. F.; Vicic, M. J. Fluid Mech. 2002, 456, 239. (6) Torquato, S. Random Heterogeneous Materials; Springer: New York, 2002. (7) Dobry, A. J. Chim. Phys. 1953, 50, 507. (8) Simha, R. J. Appl. Phys. 1952, 23, 1020. (9) Happel, J. J. Appl. Phys. 1957, 28, 1288. (10) Sherwood, J. D. Chem. Eng. Sci. 2006, 61, 6727. (11) Ruiz-Reina, E.; Carrique, F.; Rubio-Herna´ndez, F.-J.; Go´mezMerino, A.-I.; Garcı´a-Sanchez, P. J. Phys. Chem. B 2003, 107, 9528. (12) Levine, S.; Neale, G. H. J. Colloid Interface Sci. 1974, 47, 520. (13) Zharkikh, N. I.; Shilov, V. N. Colloid J. USSR 1982, 43, 865. (14) Kozak, M. W.; Davis, E. J. J. Colloid Interface Sci. 1989, 127, 497. (15) Kozak, M. W.; Davis, E. J. J. Colloid Interface Sci. 1989, 129, 166. (16) Lee, E.; Chu, J.-W.; Hsu, J.-P. J. Colloid Interface Sci. 1999, 209, 240. (17) Ding, J. M.; Keh, H. J. J. Colloid Interface Sci. 2001, 236, 180. (18) Carrique, F.; Arroyo, F. J.; Jimenez, M. L.; Delgado, A. V. J. Phys. Chem. B 2003, 107, 3199. (19) Ahualli, S.; Delgado, A. V.; Miklavcic, S. J.; White, L. R. Langmuir 2006, 22, 7041. (20) Happel, J. AIChE J. 1958, 4, 197. (21) Kuwabara, S. J. Phys. Soc. Jpn 1959, 14, 527. (22) Dukhin, S. S.; Shilov, V. N. Dielectric phenomena and the double layer in disperse systems and polyelectrolytes; Wiley: New York, 1974. (23) Fixman, M. J. Chem. Phys. 1980, 72, 5177. (24) Hinch, E. J.; Sherwood, J. D. J. Fluid Mech. 1983, 132, 337. (25) Sherwood, J. D. J. Fluid Mech. 1980, 101, 609. (26) Rubio-Herna´ndez, F.-J.; Ruiz-Reina, E.; Go´mez-Merino, A.-I. J. Colloid Interface Sci. 1998, 206, 334. (27) Sherwood, J. D.; Rubio-Herna´ndez, F. J.; Ruiz-Reina, E. J. Colloid Interface Sci. 2000, 228, 7. (28) Carrique, F.; Garcı´a-Sanchez, P.; Ruiz-Reina, E. J. Phys. Chem. B 2005, 109, 24369. (29) Ruiz-Reina, E.; Garcı´a-Sanchez, P.; Carrique, F. J. Phys. Chem. B 2005, 109, 5289. (30) Batchelor, G. K. J. Fluid Mech. 1970, 41, 545. (31) Zholkovskiy, E. K.; Adeyinka, O. B.; Masliyah, J. H. J. Phys. Chem. B 2006, 110, 19726. (32) Watterson, I. G.; White, L. R. J. Chem. Soc., Faraday Trans. 2 1981, 77, 1115. (33) Booth, F. Proc. R. Soc. London, Ser. A 1950, 203 533. (34) de Kruif, C. G.; van Iersel, E. M. F.; Vrij, A.; Russel, W. B. J. Chem. Phys. 1985, 83, 4717.
