The effect of surface charge distribution on hindered diffusion in pores

Journal of Colloid and Interface Science 338 (2009) 250–260

Contents lists available at ScienceDirect

Journal of Colloid and Interface Science www.elsevier.com/locate/jcis

The effect of surface charge distribution on hindered diffusion in pores Ronald J. Phillips * Department of Chemical Engineering and Materials Science, University of California at Davis, Davis, CA 95616, USA

a r t i c l e

i n f o

Article history: Received 3 April 2009 Accepted 27 May 2009 Available online 3 June 2009 Keywords: Hindered diffusion Electrostatic interaction Polymer gel

a b s t r a c t An analytic solution is derived for the hindered diffusion of charged, small solutes in charged, cylindrical pores in which the pore wall potential consists of the sum of an average and an oscillatory component. When the oscillatory contribution is absent, the effect of electrostatic interactions on diffusion is negligible. However, when the wall potential or surface charge density varies axially, electrostatic interactions hinder the rate of diffusion significantly, and can stop it completely if ‘‘choke points” develop where the solute concentration becomes zero. The degree of hindrance is generally weaker when the electrostatic charge on the pore wall and the charge on the solute have the same signs, leading to a repulsion, than it is in the presence of an attraction. The electrostatic hindrance is also affected by the length scale of the axial variation along the pore wall, becoming stronger as that length grows, until an asymptotic value is reached. The theory for the effect of variations of the electrostatic potential on rates of diffusion is shown to be in good agreement with experimental data taken from the literature. The results here are obtained by using generalized Taylor dispersion theory, and are therefore rigorous predictions of what occurs over times long enough that the solute diffuses through a tube many times longer than a single periodic cell. The electrostatic interactions are calculated using the linear Poisson–Boltzmann equation. Ó 2009 Elsevier Inc. All rights reserved.

1. Introduction Transport of Brownian solutes through microporous media is affected significantly by interactions between the solute and the pore surfaces. Different rates of transport caused by these interactions form the basis for many bioseparation processes, and also play a role in a number of physiological systems, and hence the subject of hindered transport has been studied for some time [1,2]. Particularly in the case of diffusion and convection of globular solutes in cylindrical pores, reliable theoretical predictions are available for conditions in which charge effects are negligible. However, when electrostatic interactions play a significant role, and when the microstructure is more complicated than that of a spherical solute in a straight, cylindrical pore, there are still a number of unanswered questions. In this paper, the problem of a charged solute in a cylindrical pore with a periodically varying wall potential is considered. The problem of how hindered diffusion in a pore is affected by a constant charge density on the wall has been treated quite recently [3]. The analysis shows that, for a spherical solute diffusing through a pore, electrostatic repulsion between the solute and the pore wall has an effect that is minor relative to the accompanying change in the equilibrium partition coefficient. Essentially, once a solute is inside a pore, the repulsive interaction causes the solute to spend more time in positions far from the wall, but * Fax: +1 916 752 1031. E-mail address: [email protected] 0021-9797/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.jcis.2009.05.063

the resulting decrease in the hydrodynamic drag is a small effect compared to the reduction in the likelihood that the solute enters the pore in the first place. As pointed out by Dechadilok and Deen [3], models of hindered diffusion in pores that assume the diffusion coefficient in the pore is equal to the value at the centerline, and hence independent of radial position, only differ from theories that account for the radial variation by approximately 10%. A simple example demonstrates that this conclusion does not extend to cases where the potential or charge density on the pore wall varies spatially. If the variation in the wall potential is great enough, there are certain to be regions of the pore that are completely inaccessible to the solute, and all macroscopic transport must cease. Depending on the spacing between these regions where the concentration approaches zero, the equilibrium partition coefficient could remain finite, and even be close to unity. Hence, variations in the wall potential must become important under some conditions, and in extreme cases yield a macroscopic diffusion coefficient of zero, even in pores where the equilibrium partition coefficient is finite and where the solute is free to diffuse locally over limited distances. Theories of hindered transport in cylindrical pores often form the basis of models of hindered transport in more complicated media, such as hydrogels, and it is noteworthy that the problem of a protein moving through a pore with varying wall potential captures some qualitative effects that are missed by the case where the wall potential is constant. In a hydrogel with a fibrous microstructure, the solute would experience gradients in the electrostatic potential in the direction of the concentration gradient,

R.J. Phillips / Journal of Colloid and Interface Science 338 (2009) 250–260

and there is no comparable electric field in a cylindrical pore with constant wall properties. Intuitively, the effect of these gradients in electrostatic potential would be alternately to speed up, and then slow down, the transport of the solute. However, as is shown here, and can be seen by the example just discussed, the symmetry of the electric field does not imply that it has no net effect on the macroscopic diffusion coefficient. The framework of generalized Taylor dispersion theory provides a rigorous and convenient method for evaluating the effect of a spatially varying potential on a pore wall. The term ‘‘Taylor dispersion” or ‘‘Taylor–Aris dispersion” is commonly used to refer to the effect of the interaction between molecular diffusion and the velocity field in a cylindrical tube on the macroscopic, or ‘‘longtime,” dispersion of the solute. Taylor [4] showed that local diffusion between positions in the fluid with different velocities causes the solute to disperse at a rate that depends on the Peclet number Pe, where

Pe ¼

u0 R : D

ð1Þ

For Poiseuille flow in a tube with maximum velocity u0, the dispersion is characterized by a dispersion coefficient Ddis that is related to the molecular diffusion coefficient D by

Ddis Pe2 ¼1þ : D 192

ð2Þ

This result agrees very well with experimental observations, and is even used in some instruments for measuring diffusion coefficients [5]. Aris [6], Brenner [7] and Brenner and Gaydos [8] re-derived and extended Taylor’s [4] result with more rigorous methodologies, and Brenner and Adler [9] generalized it further still to be applicable to any spatially periodic porous medium. Although Brenner and Gaydos [8] explicitly allow for the possibility of a constant electrostatic potential on the pore wall, the effects of a periodic variation in that potential are more easily taken into account by using the work of Brenner and Adler [9]. The theory of Brenner and Adler [9] has already been applied to the problem of hindered transport in a fibrous medium in which hydrodynamic interactions between the solute and the fibers lead to a spatially varying local diffusion coefficient [10,11]. Those results, and the underlying theory as discussed in the next section, provided motivation for describing rates of hindered diffusion in fibrous media in terms of a ‘‘steric factor” that accounts for the irregular path of the solute, and a ‘‘hydrodynamic factor” that accounts for the enhanced drag experienced in the fibrous porous medium [12,13]. In the absence of charge effects, this approach yields predictions in good agreement with experimental data for a wide array of solutes and gels [14]. Although the current work does not directly involve fibrous media, it is motivated in part by experimental results that have been reported for the diffusion of charged solutes in charged agarose gels [15,16]. Johnson et al. [15] measured rates of diffusion of negatively charged bovine serum albumin, ovalbumin and lactalbumin in agarose gel that had been modified to contain sulfate groups that carry a negative charge. Although they emphasize the greater relative importance of the partitioning effect, Johnson et al. [15] did find a significant effect of charge on rates of diffusion also, particularly for the lactalbumin. Fatin-Rouge et al. [16] used fluorescence correlation spectroscopy (FCS) to measure rates of diffusion in unmodified, and consequently weakly negatively charged, agarose. They found rates of diffusion of positively charged solutes that depended strongly on ionic strength, even in the absence of binding of the solute to the gel fibers. In the sections that follow, the mathematical formalism whereby generalized Taylor dispersion theory allows the long-time, macroscopic diffusion coefficient to be calculated from local or

251

short-time diffusion coefficients is briefly reviewed. The equations relevant to the problem of diffusion in a cylindrical pore with a periodically varying wall potential are then defined. In most instances, physical insight as well as quantitatively accurate results can be obtained most simply by radially averaging the governing equations, a process that involves approximations that are justified by comparison with more complete, numerical calculations for some test cases. The results are shown to be in good quantitative agreement with experimental data reported by Fatin-Rouge et al. [16]. 2. Generalized Taylor dispersion theory We begin with a short synopsis of the method for calculating macroscopic diffusion coefficients, known as generalized Taylor dispersion theory [8,7,9], simplified to suit the current application. Consider a charged, spherical solute with radius a diffusing in a cylindrical pore with radius R, as shown in Fig. 1. The electrostatic potential on the pore wall has both an average component w0 as well as a periodic, oscillatory component w0 . The wavelength of the oscillation in w0 is L. The particle is released at some position in the pore, and the flux vector J of the probability density P(x,t) of finding the particle at position x at time t is of the form

J ¼ UP  D  rP;

ð3Þ

where U is the solute velocity and D the diffusivity tensor. If it is transport over macroscopic length scales, much larger than L, that is of primary concern, then it may be desirable to find appropriate averages of U and D that describe so-called ‘‘macrotransport” over macroscopic distances. These averages are the mean, long-time dispersivity tensor D* and the mean, long-time convective flux U*. The averaged properties are referred to as ‘‘long-time” dispersion and convection coefficients because they apply only in the limit

L2 t >>   ; D

ð4Þ

which has the physical significance that the solute must have had ample time to sample all accessible positions in the pore. Here the term ‘‘average” does not have the same meaning as in the context of statistical thermodynamics, because the system is not at equilibrium. Simply weighting the diffusivity by the probability density and integrating, for example, would yield a finite long-time dispersion coefficient even if portions of the tube were rendered inaccessible by the existence of a high electrostatic potential. In order to compute the long-time transport coefficients, Brenner and co-workers [7,9] define a function that is the probability of finding the solute at a particular local position r in any of the unit cells, where a unit cell in this case is the section of pore in Fig. 1 in

Fig. 1. Cylindrical pore with radius R and wall potential w0 + w0 that varies with period L.

252

R.J. Phillips / Journal of Colloid and Interface Science 338 (2009) 250–260

^ It the region 0 6 z 6 L. Here we denote this probability function P. can be found by summing P over all periodic cells in the tube. The corresponding flux vector is ^J, where

^J ¼ U P ^ ^  D  rP;

ð8Þ

where B is spatially periodic in a spatially periodic medium. If no obstructions are present, then B is a linear function of position, and the concentration change C0 is similarly a linear function of position over a periodic cell. However, when obstructions are present, whether they be physical or imposed by an electric field, then B is relatively constant in regions where the concentration profile is approximately flat, and B changes abruptly in locations with steeper concentration gradients. ^ is obtained from Eqs. (6)–(10) Once the probability density P with ^J as defined in Eq. (5), the macroscopic convective flux vector U* can be obtained from Eq. (13). If the B field is then found by solving Eqs. (12)–(17), Brenner and Adler [9] show that the macroscopic transport coefficient D* is given by

ð9Þ

D ¼

ð5Þ

^ is governed by the equation and the probability density function P

r  ^J ¼ 0;

ð6Þ

with the boundary conditions

er  ^J ¼ 0 at r ¼ R; ^ @P ¼ 0 at r ¼ 0; @r ^ ; ^ Pjz¼0 ¼ Pj z¼L

ð7Þ

Z

^ ðrBÞt  D  rBdV: P

ð19Þ

VP

and

^ ^ rPj z¼0 ¼ rPjz¼L :

ð10Þ

Eq. (7) is a requirement that there be no flux through the pore wall, where the unit normal vector is er, while Eq. (8) reflects the fact that the centerline of the pore is a line of symmetry. The two conditions ^ and its gradient to be periodic. Finally, the (9) and (10) require P probability density is normalized,

Z

^ PdV ¼ 1;

ð11Þ

VP

which eliminates the possibility of a trivial solution to Eqs. (3), (6)– (9) and (10). The term VP in Eq. (11) is the volume of one periodic cell of the pore. ^ alone is not sufficient to enable The probability density P or P one to calculate the coefficients that describe convection and dispersion over macroscopic regions. In generalized Taylor dispersion theory, the additional information is contained in a vector field B (which has no relation to magnetism). The B-field is the solution to





^  rB  J  rB ¼ U  ; r  PD

ð12Þ

where U* is the macroscopic convective flux vector defined by

U ¼

Z

^JdV:

ð13Þ

VP

The boundary conditions for B in the current problem are

er  D  rB ¼ 0 at r ¼ R; @B ¼ 0 at r ¼ R; @r   Bjz¼L  Bjz¼0 ¼  rjz¼L  rjz¼0 ;

ð14Þ ð15Þ ð16Þ

and

rBjz¼L ¼ rBjz¼0 :

ð17Þ

The vector r in Eq. (16) denotes geometrically equivalent positions at the same (r, h) positions of the planes bounding the periodic cell at z = 0 and z = L. Because Eqs. (12)–(17) only involve gradients of B, or in the case of Eq. (16) the jump in B from one boundary to the next, B is only defined to within an arbitrary constant. Koch et al. [17] point out that the B-field can be interpreted as the vector field that relates local fluctuations in the solute concentration (or probability density) to the average of the gradient in the concentration. In other words, if a constant, average concentration gradient rhCi is imposed over a macroscopic region of a porous or fibrous medium, the concentration fluctuations C0 may be expressed in the form

C 0 ¼ B  rhCi;

ð18Þ

The superscript t in Eq. (19) indicates the transpose of the tensor in parentheses. The macroscopic coefficient U* in Eq. (13) describes how the mean position of a tracer particle introduced into the porous medium varies with time, over macroscopic distances. Similarly, the dispersion tensor D* describes how the mean-squared displacement of that tracer particle varies with time, again over macroscopic distances, and accounting for the effects of the spatially varying velocity field. In transport through a cylindrical pore over macroscopic distances, it is generally the zz component of D* that is of interest. Taking the double dot product of both sides of Eq. (19) with ezez, where ez is a unit vector in the z direction, and making use of the simplification

D ¼ DI

ð20Þ

that is valid for small solutes for which a/R  1, yields the result

Dzz ¼ D

Z

^ rBz  rBz dV: P

ð21Þ

VP

The coefficient Dzz in Eq. (21) describes the growth in the meansquared displacement of a tracer in a pore, in the presence of effects due to velocity or electric fields. In the sections below, Dzz is referred to as a macroscopic diffusion coefficient, rather than dispersion coefficient, because the term ‘‘dispersion” is so often reserved for situations involving fluid flow. ^ For spherical Brownian solutes, the probability density PðrÞ refers to the probability of locating the center of the solute at a particular position r within one of the periodic cells of the pore. The velocity U is the local solute velocity, and is in general not equal to the local fluid velocity at r. One reason why the particle velocity can differ from the fluid velocity is hydrodynamic interactions between it and the pore wall. In addition, if electrostatic interactions exert a force F, the result is a solute velocity of the form MF that must be included in U, where M is the solute mobility that is related to the local diffusivity tensor by the Einstein relation

M¼

D : kT

ð22Þ

Here k is the Boltzmann constant and T is temperature. If hydrodynamic interactions are negligible, then the mobility tensor for a spherical solute with radius a in a liquid with viscosity l is given by Stokes drag law,

M ¼ nI;

ð23Þ

where n = 1/(6pla). Provided the electrostatic force is accounted for in computing the solute velocity U, the theoretical results above may be applied rigorously to analyze the transport of charged solutes through pores with charged walls [9].

R.J. Phillips / Journal of Colloid and Interface Science 338 (2009) 250–260

In this development, the approximation is made that the solute velocity U is given by

U ¼ qenrw;

ð24Þ

where q is the charge on the solute, e the unit of electrostatic charge, and w the electrostatic potential. Eq. (24) is appropriate when the solute is small relative to the Debye length [18], and the surface interaction between the solute and the pore wall is negligible (i.e., the surface potential of the solute does not change in response to the field created by the pore wall, and vice versa). It has also been assumed that the particle radius is small enough relative to the pore radius that hydrodynamic interactions with the wall are negligible. The analysis here therefore isolates electrostatic effects relative to hydrodynamic effects, in a straight pore where the tortuosity plays no role. 2.1. Relation to steric and hydrodynamic factors in fibrous gels For hindered diffusion in uncharged, cylindrical pores with diameter comparable to that of the diffusing solute, it is mentioned above that the hydrodynamic effect can be accounted for with an error of less than 10% by using hydrodynamic mobilities calculated along the pore centerline. In other words, the mobility of the solute may be affected significantly by being located in a pore with comparable dimensions, but the local position within the pore is not of great significance. The same is generally true in fibrous media: the fact that solute motion requires displacement of fluid through an effective porous medium is more important than the specific position relative to the nearest fiber. Given the geometric complexity of many fibrous media, and the many approximations that must be made to model hindered transport in general, it is therefore reasonable to approximate the local diffusion tensor D in Eq. (19) with a spatially averaged result that accounts for hydrodynamic interactions only in an averaged sense. This approximation simplifies Eq. (19) to

D ¼ hDi

Z

^ ðrBÞt  rBdV: P

ð25Þ

VP

The integral term in Eq. (25) accounts for steric or tortuosity effects, and also for the interaction between diffusion and the velocity and electric fields that results in dispersion. The hydrodynamic effect is accounted for only by changing D to hDi to reflect the fact that the solute is not moving through an unbounded fluid. This separation of steric and hydrodynamic effects in modeling hindered transport has been applied successfully in uncharged systems [14]. The steric terms can be found from macrotransport theories, including Eqs. (5)–(17) above as well as other approaches [19,20], or from stochastic simulations [21]. The hydrodynamic effect can be found from the Brinkman effective medium approximation [12] or from direct calculations of solute mobilities in three-dimensional fiber arrays [13,14].

253

interest here is the case when the potential on the pore surface varies about a mean value aw0 according to

 

2p z w ¼ w0 a þ cos at r ¼ R; L

ð27Þ

and the potential is symmetric about the centerline,

@w ¼ 0 at r ¼ 0: @r

ð28Þ

In accordance with the discussion in the last section, we are concerned with ionic strengths such that the Debye lengths 1/j satisfy ja  1, and with solute radii a that satisfy a/R  1. The solution to Eq. (26) with the boundary conditions (27) and (28) is

 

I0 ðjrÞ 2pz I0 ðbrÞ ; wðr; zÞ ¼ w0 a þ cos I0 ðjRÞ L I0 ðbRÞ

ð29Þ

where I0(x) is a modified Bessel function of the first kind, of order zero, and the parameter b is given by

b¼

"  #1=2 2 2p þ j2 : L

ð30Þ

As can be seen from Eq. (29), the magnitude of the potential is a maximum at the entrance, z = 0, or exit, z = L, of a periodic cell in the pore. It is minimum at the midpoint, z = L/2. Although Eq. (29) is the solution when the surface potential is specified according to Eq. (27), it can be related to the case where the surface charge density r is specified by using



@w r ¼ n  rwjr¼R ¼   : @r

ð31Þ

r¼R

Eq. (31) shows that, like the surface potential, the surface charge density is comprised of mean and oscillatory components. The solution Eq. (29) is plotted in Figs. 2 and 3 for a pore with length L/R = 4 and a = 1. With a = 1, there is a mean surface potential equal to w0, and the oscillation from the cosine term has an amplitude equal to the mean. From Fig. 2, one can see that for a Debye length such that jR = 1, there is less variation in w across the pore than might be expected from the exponential decay that occurs near a flat surface. The maximum variation occurs at ^z ¼ z=L ¼ 0, where the contribution of the cosine term is greatest, but even there the change is less than a factor of 2. The change along the length of the pore is shown in Fig. 3, for three radial positions. The three curves are qualitatively similar, showing maxima at z/R = 0 and z/R = 4 and decaying to zero at the midpoint, z/R = 2. Because Eq. (26) is linear, the sum of any number of solutions such as those shown in Figs. 2 and 3 is also a solution, by

3. Effect of uneven pore wall potential If a pore is filled with a 1–1 electrolyte, and the potential on the pore walls and solute surface are comparable to or less than kT/e, where e is the unit of electrostatic charge, then the electrostatic potential in the liquid-filled pore should be in good qualitative, and in most cases quantitative, agreement with the prediction of the linear Poisson–Boltzmann equation [22],

r2 w ¼ j2 w;

ð26Þ

where j is the inverse of the Debye length [23]. (For a 1–1 electrolyte, the linearization used to obtain Eq. (26) amounts to setting sinh(x)  x, an approximation that is in error by 18% at x = 1.) Of

Fig. 2. Electrostatic potential w vs. radial position r/R at three axial positions, ^z ¼ z=R, of a periodic cell with L/R = 4, jR = 1 and a = 1.

254

R.J. Phillips / Journal of Colloid and Interface Science 338 (2009) 250–260





^  rB ¼ 0: r  PD

ð37Þ

Expanding the r operators in Eq. (37) in terms of the cylindrical coordinate system shown in Fig. 1, making use of Eq. (20), and dotting each term in the equation by the constant unit vector ez, then yields

   

@ ^ @Bz 1 @ ^ @Bz P D rP þ ¼ 0: @z r @r @z @r

ð38Þ

Eq. (38) is to be solved subject to the boundary conditions (14)– (17), which simplify to

Fig. 3. Electrostatic potential w vs. axial position z/R at three radial positions, ^r ¼ r=R, of a periodic cell with L/R = 4, jR = 1 and a = 1.

the principle of linear superposition. In this way, a Fourier cosine representation of an arbitrarily complicated wall potential could be formed, and the methods used here to obtain Dzz would still be applicable. 3.1. The probability density Because the electrostatic contribution to the solute velocity is proportional to rw as shown in Eq. (24), in the absence of flow a ^ can be found quite readily. solution for the probability density P By inspection one can see that ^J ¼ 0 is a solution to Eqs. (5)–(10). Then a possible solution may be found by equating the convective and diffusive contributions to ^J,

^ rw ¼ DrP: ^ qenP

ð32Þ

In light of the Einstein equation (22), which simplifies to D = kTn if hydrodynamic interactions are negligible, the solution to Eq. (32) is

^ ¼ K p eqew=kT ¼ K p eE^w^ : P

ð33Þ

^ is the dimensionless electrostatic potential Here the parameter E energy of the solute,

^ ¼ qew0 ; E kT

ð34Þ

^ is the dimensionless electrostatic potential, w

^¼ w w w0

ð35Þ

and Kp is

Kp ¼

Z

qew=kT

e

1 dV

;

ð36Þ

VP

@Bz ¼ 0 at r ¼ R; @r @Bz ¼ 0 at r ¼ 0; @r Bz jz¼L  Bz jz¼0 ¼ L;

ð39Þ ð40Þ ð41Þ

and

rBz jz¼L ¼ rBz jz¼0 :

ð42Þ

Once obtained, the solution to Eqs. (38)–(42) can be used in conjunction with Eq. (21) to obtain the macroscopic diffusion coefficient Dzz . The fact that the radial variations in Figs. 2 and 3 are modest compared to the axial variations suggests that a useful simplification can be achieved by averaging over the radial coordinate. Let the radial average hSi of S be defined by

hSi ¼

2 R

2

Z

R

Srdr;

ð43Þ

0

^ or some combination of the two. The averwhere S could be Bz ; P, aged probability density function, found by integrating the solution in Eq. (33) using the potential from Eq. (29), is plotted in Fig. 4. Averaging over the radial coordinate in this fashion permits an analytical solution to be obtained, without sacrificing the important physical characteristics of the problem. The quantitative effect of pursuing a solution via this averaging method is assessed by direct numerical evaluation in below. Averaging the terms in Eq. (38) in accordance with Eq. (43), one obtains

   Z R d ^ @Bz 2 @ ^ @Bz dr ¼ 0: P þ 2 rP dz @z @r R 0 @r

ð44Þ

The integral of the r-derivative term in Eq. (44) is trivial, and use of Eqs. (39) and (40) eliminates that term. A second integration of the remaining term gives

 ^ @Bz ¼ K 2 ; P @z

ð45Þ

as required by the normalization condition (11). The dimensional potenial w is given by Eq. (29). In the absence of flow, Eq. (33) shows that the probability of finding the solute at a particular position in a periodic cell is governed by a Boltzmann distribution. In ^ shows the relight of Eq. (33), the normalized probability function P verse of the trends of the potential w. 3.2. The B field To derive the solution for B in the absence of flow, we first simplify Eq. (12) by noting that, as just shown, ^J ¼ 0 if the fluid velocity is zero. The convective term proportional to rB in Eq. (12) may therefore be eliminated. In addition, according to Eq. (13), the macroscopic convective flux U* is found by integrating the flux ^J over the volume of the pore VP, and is also zero. Eq. (12) therefore immediately simplifies to

^ vs. axial position z/R for a periodic Fig. 4. Radially averaged probability density hPi cell with L/R = 4 and a = 1.

R.J. Phillips / Journal of Colloid and Interface Science 338 (2009) 250–260

where K2 is a constant. Making the approximation

 @Bz ^ dhBz i ; P^  hPi dz @z

ð46Þ

one finds that the Bz derivative is inversely proportional to the average probability density,

dhBz i K 2 ; ¼ ^ dz hPi

ð47Þ

^ in Eq. (47) is known from Eq. The average probability density hPi (33), and the solution (47) satisfies the periodicity requirement (42) because of the periodic nature of w. The constant K2 is specified by Eq. (41), and is given by

K2 ¼ R L

L

^ 1 dz hPi 0

ð48Þ

;

with the magnitude of hBzi following by straightforward integration of dhBzi/dz. Since it is changes in B that have physical significance, not the magnitude itself, the value assigned to hBzi at the edge of a periodic cell is arbitrary. The complete solution for the radial average hBzi is therefore specified by Eqs. (33), (36), and (47) and (48). The derivative is of most importance, and that is given by

dhBz i ¼ h RL RR dz 0

0

eqew=kT rdr

L i1 nR o; R qew=kT dz e rdr 0

ð49Þ

where w is specified by Eqs. (29) and (30). The z-dependence of dhBzi/dz determines the effect of local constrictions caused by locally high values of the potential. In Eq. (49), the z-dependence derives entirely from the second term in the denominator, which is not integrated over the pore length. At positions of high potential, the exponential becomes very small, leading to locally steep gradients in hBzi(z). ^ and dhBzi/dz has signifiThe inverse relationship between hPi cant implications for the macroscopic dispersion coefficient Dzz , which in terms of the radially averaged quantities takes the form (cf. Eq. (21))

hDzz i ¼ pR2 D

Z 0

L

^ hPi

 2 dhBz i dz: dz

ð50Þ

From Eq. (47), in any pore with positions where the probability density becomes very low, the z-deriviative of Bz must become very high. However, the overall change in Bz from one end of a periodic cell to the other is specified by the boundary condition (41). The steep gradients must therefore be constrained to a very small ^ and region of the pore, and from Eq. (50), it is the solution for hPi hBzi in that region that largely determines the macroscopic coefficient Dzz . As explained in Section 2, the B field relates local concentration (or probability density) fluctuations in a pore or porous medium to the overall macroscopic concentration gradient that is causing the flux. The presence of an abrupt change in B therefore indicates an abrupt concentration fluctuation, which in the current example occurs at any location where the local probability density becomes very small. For a pore in which there are regions of very high electrostatic potential separated by regions of low potential, the implication is that the solute concentration is relatively constant in the low-potential region, and goes through very sudden changes in the high-potential regions. As the primary barrier to transport, the characteristics of the high-potential regions determine the macroscopic coefficient Dzz , and these regions are therefore given extra weight in the integral specified in Eq. (50). Results for Bz in charged pores are plotted in Fig. 5 for solutes ^ ¼ 1; 2 and 3. The pore with electrostatic potential energies of E

255

length is L/R = 4, and the Debye lengths are such that jR = 1. The relatively abrupt changes in Bz that are apparent near z/R = 0 and z/R = 4 indicate significant fluctuations in the probability density in those regions away from its unperturbed, linear profile. The effect becomes significantly stronger for solutes with higher charge or potential energy, such that Bz is almost unchanging in the cen^ ¼ 3, indicating a similarly unchanging tral region of the pore for E concentration as indicated by Eq. (18). The physical interpretation of Fig. 5 is therefore that, when the potential is high at the entrance and exit of the periodic cell, the concentration there changes very rapidly. In the central region where the potential is relatively low, the concentration is relatively high and approximately constant in magnitude. There exist, in a sense, concentrated pools of solute separated by highly restrictive entrances and exits. In such a scenario, one would expect to find very different behavior on short and long time scales, because over short time scales the solute is never required to traverse the highpotential regions. Brownian dynamics simulations of solutes in charged networks do show a transition as long-time, asymptotic behavior is approached [24]. In a pore such as those considered here, anomalous diffusion would not be present at long times, as has been suggested for other systems [25]. However, it would take progressively longer and longer times for the asymptotic behavior to be observable, since many periodic cells must be sampled by the solute tracer. The true macroscopic diffusion coefficient could become virtually unattainable in dynamic simulations, if it is reduced by orders of magnitude by electrostatic hindrance. It is the derivative dhBzi/dz that determines the macroscopic dispersion coefficient, and that quantity is plotted in Fig. 6. The derivative, which is squared in the integral term of Eq. (50), is nearly zero in the central region of the periodic unit cell. The averaging is therefore heavily weighted toward the entrance and exit. If the local molecular diffusivity D were a function of axial position, and not constant, the result would be an overwhelming emphasis on its value in the region near z/R = 0. However, even in the current case where D is constant, the integral in Eq. (50) shows a significant effect from the electrostatic interaction, as discussed in next section. 3.3. Results for Dzz Values for the macroscopic, long-time dispersion coefficient Dzz are shown in Fig. 7 for a tube with periodic cells of length L/R = 4 in which a = 1. Results are presented for dimensionless solute poten^ ¼ 1; 2 and 3, for Debye lengths in the range tial energies of E 0.1 6 jR 6 7.5. There is clearly a strong effect, particularly at the ^ ¼ eqw =kT. The electrostatic effect is strongest higher values of E 0

Fig. 5. Radially averaged B field hBzi vs. axial position z/R for a periodic cell with L/ R = 4 and a = 1. Results for three values of the dimensionless solute potential energy ^ are shown. E

256

R.J. Phillips / Journal of Colloid and Interface Science 338 (2009) 250–260

Fig. 6. Derivative of the radially averaged B field dhBzi/dz vs. axial position z/R for a periodic cell with L/R = 4 and a = 1. Results for three values of the dimensionless ^ are shown. solute potential energy E

^ ¼ 1; 2 and 3 in a periodic Fig. 8. Macroscopic diffusion coefficient Dzz vs. jR for E cell with L/R = 4 and a = 0.

^ ¼ 1; 2 and 3 in a periodic Fig. 7. Macroscopic diffusion coefficient Dzz vs. jR for E cell with L/R = 4 and a = 1.

again, the results in the limit jR ? 0 are unchanged in all three Figs. 7–9. The similarities and differences in these figures may be interpreted by noting that the main influence of the presence of a constant, mean potential on the pore wall is to move the solute, on average, to positions closer or farther from it. When the solute and the wall undergo a repulsive interaction, the solute is pushed to the center of the pore. If the Debye length is comparable to or smaller than the pore radius, then electrostatic effects are significantly weaker at the pore center, and screening is therefore very effective. The only exception occurs in the limit jR ? 0, in which case electrostatic effects are as strong at the pore center as they are near the wall. When the solute–wall interaction is attractive, the effect is to move the particle closer to the wall, where electrostatic effects are strong even for small jR, hindering the effect of Debye screening. The case where the mean wall potential is zero falls between the other two. These results therefore indicate that increasing the salt concentration, which reduces the dimension of the double layer, reduces the electrostatic effect much more effectively for repulsive interactions than for attractive ones. This conclusion is clearly evident in Fig. 10, where D*/D is shown when the pore surface and solute undergo a mean repulsive interaction (a = +1), no interaction (a = 0), ^ ¼ 2. Even when the Deand an attractive interaction (a = 1), for E bye length is small, jR = 10, the electrostatic effect is strong when the mean solute–wall interaction is attractive. Brownian dynamics simulations have also shown an enhanced electrostatic hindrance in fibrous gels when the solute and gel fibers have opposite charges, relative to the similar charge case, in the limit where the gel fibers are dilute [24,26]. The limit of dilute fiber concentration is comparable to the limit a/R  1 considered here.

when the Debye length is very large relative to the tube radius, so that jR ? 0, but even at jR = 5 it reduces the rate of axial disper^ ¼ 3. The effect becomes small sion by approximately 20% when E when jR > 7.5. The effect of using radial averaging to simplify the solution to Eq. (37) was checked by obtaining the complete solution numerically, and comparing the two sets of results for Dzz . To obtain the numerical solution, the derivatives in Eq. (38) were approximated by central difference formulas, and the equation was then solved for Bz iteratively by using a Gauss–Seidel method [27], using ‘‘fictitious” nodes for derivative conditions at the boundaries of the peri^ ¼ 1 and E ^ ¼ 3, at Debye odic cell [28]. Calculations were done for E lengths corresponding to jR = 1, 2 and 4. The results, shown as the solid points in Fig. 7, indicate that the radial averaging has a negligible impact on the calculation of Dzz . It therefore appears that Eqs. (49) and (50) give an excellent approximation to the full numerical solution. The electrostatic hindrance persists, and even grows stronger, when the mean potential on the pore wall is zero, a = 0, but the oscillatory contribution remains. Results for that case are shown in Fig. 8, for conditions that are otherwise the same as in Fig. 7, and show very similar qualitative behavior. In fact, in the limit of very large Debye lengths, jR ? 0, the results in Figs. 7 and 8 are the same. However, the screening effect that is present at finite salt concentrations is relatively stronger when the interaction between the pore wall and solute is repulsive. The persistence of the electrostatic hindrance is magnified further still when there is an attraction between the solute and the pore wall, as shown in Fig. 9. There, the rate of axial diffusion is reduced by approximately 20% ^ ¼ 2, and by more than 50% for E ^ ¼ 3. Once even at jR = 15 for E

^ ¼ 1; 2 and 3 in a periodic Fig. 9. Macroscopic diffusion coefficient Dzz vs. jR for E cell with L/R = 4 and a = 1.

R.J. Phillips / Journal of Colloid and Interface Science 338 (2009) 250–260

257

^ ¼ 2 in a Fig. 10. Macroscopic diffusion coefficient Dzz vs. jR for a = 1, 2 and 3 with E periodic cell with L/R = 4.

Fig. 12. Macroscopic diffusion coefficient Dzz vs. the dimensionless periodic cell ^ ¼ 1; 2 and 3 and a = 1. length L/R for E

If the surface charge density r is held constant while jR is varied, then the surface potential must vary in accordance with Eq. ^ must vary as specified by Eq. (31), and in turn the parameter E (34). For a positively charged surface, if jR decreases and r is con^ increase. The effect of holding stant, the surface potential w0 and E ^ ¼ 1 and E ^¼3 the surface charge density fixed on the curves for E in Fig. 7 is shown in Fig. 11, for conditions where a = +1 and L/R = 4. The decrease in the surface potential as jR increases has a strong effect, causing an abrupt rise in D*/D. The effect of electrostatic interactions is reduced to being negligible at jR  4, even for ^ ¼ 3. In addition, the steepest changes in D*/D occur at the smallE est values of jR, in contrast with the curves in Fig. 7, where D*/D approaches a clear limit as jR ? 0. Note that, because the surface interaction between the solute and pore wall is neglected, any of the results in Figs. 7–10 pertain to either constant-charge-density or constant-potential boundary conditions on the pore wall. It is only the question of what is held constant, potential or charge density, as jR changes, that differs between those figures and Fig. 11. The length of the periodic cell, which corresponds to the wavelength in the cosine term of Eq. (27), also strongly effects the strength of the effect of an uneven wall potential. As shown in Fig. 12, in the limit of small wavelength, where L/R ? 0, the electrostatic effect becomes negligible, and the solution reduces to ^ The hindrance effect contribDzz ¼ D regardless of the value of E. uted by the changing wall potential increases as the length of the cell L/R increases, and asymptotically approaches a limiting value for L/R > 10. This behavior may be explained by noting that the electrostatic hindrance is essentially a convective effect, in which a solute velocity is driven by the electrostatic force, and is represented by

first-order derivatives in Eqs. (6) and (12). By contrast, diffusive transport is represented mathematically by second-order spatial derivatives. The physical consequence of this mathematical difference is that diffusion is relatively stronger over shorter length scales, and convection relatively stronger over larger length scales. The fluc^ that occur at very small values of L/R are relatively weak, tuations in P ^ is because diffusion is strong enough to dissipate them. When P approximately constant, Bz varies linearly, and macroscopic dispersion is governed by the molecular diffusion coefficient D. However, diffusion is a much weaker effect when L/R = 6. Therefore when w ^ are not diminished changes over longer distances, variations in P by diffusion, and variations in Bz and Dzz result. 3.4. Comparison with experimental data Although there are relatively few experimental studies of hindered diffusion designed specifically to measure electrostatic effects, some of the experiments of Fatin-Rouge et al. [16] are well-suited for comparison with the current theory. Fatin-Rouge et al. [16] used fluorescence correlation spectroscopy (FCS) to measure rates of diffusion of rhodamine 6G (R6G2+) in 1.35% agarose gels at different concentrations of NaCl. Although they also report data for other solutes, in their paper they single out the data for R6G2+ as being particularly reliable, because with R6G2+ they were able to measure rates of diffusion directly in the gel by FCS, without the need for any correction factors. The small size of R6G2+, relatively large pore sizes of agarose gels, and weak charges on the surface of agarose fibers make their data particularly relevant to this work. Fatin-Rouge et al. [16] give the hydrodynamic radius of R6G2+ as 0.744 nm, while neutron scattering experiments reveal pore radii in their agarose gels to be in the range 35–47 nm. Based on the physical structure and chemical composition of the agarose fibers, they estimate the area per surface charge to be 24.4 nm2, with the charges distributed evenly over the fiber surface. This area per surface charge corresponds to a surface charge density of  0.00656 C/m2. The solution to Eq. (26) for the potential around a charged fiber with radius Rf in an electrolyte solution is

wðrÞ ¼ wf

K 0 ðjrÞ ; K 0 ðjRf Þ

ð51Þ

where wf is the potential on the fiber surface and K0 is a modified Bessel function of the second kind, of order zero. From Eqs. (51) and (31), allowing for the change in sign of the normal vector, the surface charge density rf on the fiber is given by Fig. 11. Effect of constant-charge-density boundary conditions on the macroscopic ^ ¼ 1 and diffusion coefficient Dzz . Curves (b) and (d) are identical to the curves for E ^ ¼ 3 in Fig. 7. E

rf ¼ wf j

K 1 ðjrÞ : K 0 ðjRf Þ

ð52Þ

258

R.J. Phillips / Journal of Colloid and Interface Science 338 (2009) 250–260

To obtain parameters for use with the pore theory above, the surface charge density of 0.00656 C/m2 was substituted into Eq. (52), and the surface potential wf was calculated assuming the fiber radius Rf is 1.5 nm, as specified by Fatin-Rouge et al. [16]. The resulting surface potential is 15.7 mV at an ionic strength of 0.01 M, and is 28.8 mV at an ionic strength of 0.001 M. These potentials are low enough that the linearized Poisson–Boltzmann equation is acceptable, given the other approximations that are unavoidable in comparisons of this type. (Note, for example, that a fiber radius of 1.9 nm has also been given for agarose gels [15].) In order to predict the electrostatic effect on the rate of diffusion by using Eq. (50), it is necessary to estimate the parameters Rp, j, and L, in addition to the surface potential w0 or surface charge density r. We used a pore radius of Rp = 40 nm as being in the middle of the range given by Fatin-Rouge et al. [16]. The Debye length was calculated from its definition [23]

j¼

 2 1=2 2e I ; kT

of 0.9 can be obtained by observing the limiting value of the experimental data in conditions where electrostatic effects are completely screened, in this case corresponding to I > 0.01 M. Alternatively, it can be calculated from existing theories for hindered diffusion in the absence of electrostatic effects. Phillips [14] collected theoretical results from various sources, including Johansson et al. [21], Johnson et al. [12] and Clague and Phillips [13], and proposed the model

D b 1:09 ¼ e0:84f ea/ ; D

ð54Þ

where k is the ratio of the fiber radius to the solute radius, k = Rf/a, / is the fiber volume fraction, the parameter f is defined by

f ¼



1þ

1 k

2 /;

a ¼ 3:727  2:460k þ 0:822k2 ;

ð55Þ ð56Þ

ð53Þ

where I is the ionic strength. The ionic strengths were calculated by summing the minimum value of 8.8  104 M that Fatin-Rouge et al. [16] state is contributed by the counterions of the agarose polymer and the added salt [NaCl]added given in their Fig. 8. To ob^ we used w0 = wf/2 with tain the surface potential, and hence E, a = 1, so that the charge on the pore wall oscillates between being equal to that of the agarose fibers, which bound the ‘‘pores” or constrictions in the gel, and being zero. Finally, the ratio L/R was taken to be large, L/R ? 1, because of the expectation that the significant constrictions in the gel are sufficiently far apart that concentration fluctuations are not diminished by the influence of diffusion. Results for D*/D calculated using these parameters are plotted in Fig. 13, along with the experimental data of Fatin-Rouge et al. [16]. The theoretical calculation was done for three sets of conditions. In one case, a solute with charge q = +2 was used, in keeping with the valence of +2 on the R6G2+ solute. For comparison, the same calculation was done with a solute having a valence of 2. The prediction for the latter case shows a significantly weaker effect because, as discussed in Section 3.3, Debye screening is less effective when an attractive interaction causes the solute to prefer locations close to the charged wall. Finally, a constant-potential wall boundary condition was considered, using an intermediate potential that corresponds to a surface charge density of 0.00656 C/m2 at jR = 6. The abrupt change of D*/D with jR is not captured by the solution with a constant wall potential, indicating that the constant-charge-density boundary condition is more appropriate. The effects of hydrodynamic and steric interactions, which are not included in the pore theory here, are accounted for in Fig. 13 by multiplying D*/D as calculated from Eq. (50) by 0.9. The value

and

b ¼ 0:358 þ 0:366k  0:0939k2 :

ð57Þ

The first exponential term on the right side of Eq. (54) accounts for steric effects, and is a correlation of the simulation results of Johansson et al. [21] that was proposed by Johnson et al. [12]. The second term is a correlation of calculations by Clague and Phillips [13] of the effect of a random array of fibers on the hydrodynamic drag on a spherical solute (i.e., hDi/D in Eq. (25)). If the radius of R6G2+ is taken as 0.744 nm and the fiber radius in agarose gel is 1.5 nm (i.e., k = 2.02), then at / = 0.0135 Eqs. (54)–(57) predict D*/D = 0.89, differing from the value of 0.9 reported by Fatin-Rouge et al. [16] by only about 1%, with no adjustable parameters. Interestingly, the steric term in Eq. (54) is 0.982, while the hydrodynamic term is 0.907, indicating that the hydrodynamic effect is some five times stronger than the steric effect under these conditions. As mentioned in Section 1, Johnson et al. [15] also report experimental data related to the effect of charge of diffusion through agarose gels. Their solutes were the globular proteins bovine serum albumin (BSA), ovalbumin and lactalbumin, which had net charges (including the charge on the attached fluorescein) of 37, 19 and 9, respectively. The agarose was in the form of sulfated Sepharose beads, with an estimated surface charge density of 0.42 C/m2. The diffusion coefficient of the BSA in the gel, measured by fluorescence recovery after photobleaching (FRAP), showed no effect of electrostatic charge over ionic strengths ranging from 0.01 M to 1 M. The influence of electrostatic interactions was apparent, but weaker than expected, for the other two solutes. That the effect of repulsive electrostatic interactions between similarly charged solutes and fibers would be weaker than attractive electrostatic interactions is consistent with the results presented here. However, the surface charge density in the experiments of Johnson et al. [15] was 64 times greater than that in the experiments of Fatin-Rouge et al. [16], and the charge on the solutes was also much greater. The surface potentials were therefore so great that neither the linear nor the non-linear Poisson–Boltzmann equation could be expected to be reliable guides for interpreting the data, and a molecular-level theory would most likely be needed. 3.5. Relative influence of partitioning and hindered diffusion

Fig. 13. Comparison between predicted diffusion coefficients Dzz and experimental data for R6G2+ reported by Fatin-Rouge et al. [16].

In addition to their influence on rates of diffusion, electrostatic interactions strongly affect the equilibrium distribution of charged solutes between bulk solution and a pore. For solutes small compared to the pore radius, for which steric exclusion from the space near the pore wall is not important, the equilibrium partition coefficient U is given by

R.J. Phillips / Journal of Colloid and Interface Science 338 (2009) 250–260

U¼

1 VP

Z

eqew=kT dV:

ð58Þ

wðRÞ ¼

w0

if 0 6 z 6 ‘;

0

if ‘ 6 z 6 L:

ð59Þ

The region with length ‘ over which the pore wall is charged is short compared to the uncharged length, or ‘  L. Furthermore, both sections of the pore are long compared to the Debye length 1/j, so that to a good approximation the potential w‘ in the charged region is

w‘ ðrÞ ¼ w0

I0 ðjrÞ I0 ðjRÞ

for 0 6 z 6 ‘;

ð60Þ

whereas w(r)  0 for ‘ < z < L. A radially averaged Boltzmann factor hG‘i may be defined as

2 R

2

Z

R

eqew‘ =kT rdr;

ð61Þ

hG‘ i

  2 L  1: ‘

ð62Þ

In other words, the charged region represents a short, highly charged ‘‘choke point” in an otherwise uncharged pore. Substitution of the solution in Eq. (60) into Eqs. (49) and (50) to obtain Dzz , and into Eq. (58) to obtain U, yields

hDzz i ¼ hG‘ i

  L D ‘

nqe

Dw0 kTn ; > R R

ð65Þ

or when

qeDw0 > 1: kT

ð66Þ

Even in a medium with no net charge, if the variations in the charge density are such that the criterion (66) is satisfied, electrostatic interactions can be significant. If the medium physically resembles a continuous, solid phase permeated by a network of pores, then local electrostatic effects could reduce the rate of diffusion by orders of magnitude, even if the partition coefficient remains O(1). Although the more open pore structure of many fibrous media makes it unlikely that macroscopic diffusion would be effectively stopped, it could still be slowed considerably. The analytic results presented here provide a useful tool for estimating when electrostatic hindrance might play an important role in determining rates of diffusion, and in some cases can yield quantitative predictions for the magnitude of that hindrance.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]

ð63Þ [12] [13] [14] [15]

and

U ¼ 1:

When a charged solute diffuses in a microporous medium in which the charge is distributed homogeneously, and in which the Debye lengths are such that the electrostatic potential is approximately constant, the effect of electrostatic interactions on rates of diffusion is negligible. However, in a general, heterogeneous, charged microporous medium, there exist pore-scale variations in the distribution of charge that create a locally varying electric field. When a charged solute diffuses along a path in which its motion is alternately aided and opposed by such an electric field, local concentration fluctuations result, and the rate of diffusion over macroscopic distances is lowered. This electrostatic hindrance is expected to be significant whenever the local flux induced by the electric field is comparable to the rate of diffusion. For a solute with charge q and mobility n, in a medium in which the potential varies by Dw0 over pore dimensions of order O(R), the electrostatic hindrance is significant relative to the diffusive flux when

0

and the wall potential in the region 0 6 z 6 ‘ is taken to be high enough so that



4. Conclusions

VP

If the electrostatic potential on the pore wall is constant over the entire pore length, then under the conditions considered here the long-time diffusion coefficient is equal to the molecular diffusion coefficient, Dzz ¼ D, and it is only the partitioning effect that reduces rates of transport from one end of a pore to the other. However, if the variation in the wall potential is significant, then the two effects ^ ¼ 1; 2 and 3 with a = jR = 1 can be comparable. For example, for E and L/R = 4, partition coefficients U calculated from Eq. (58) are equal to 0.47, 0.27 and 0.18, respectively. Comparing with the results for Dzz in Fig. 7, one finds that the effects of the electrostatic interactions on diffusion and partitioning are both important with repulsive interactions. If the interaction is attractive, a = 1, but the pore is otherwise unchanged, then the partition coefficients ^ ¼ 1; 2 and 3 are 2.80, 10.1 and 43.2, respectively. With an for E attractive interaction the partition coefficient therefore increases by more than the diffusion coefficient decreases in Fig. 9, but the difference is less than an order of magnitude. As stated in Section 1, for some types of pores it is possible for ‘‘choke points” to be present such that the rate of diffusion is reduced almost to zero, while the partition coefficient remains close to unity. As an example, consider a pore with radius R and periodic cells of length L in which the wall potential is given by

hG‘ i ¼

259

ð64Þ

To obtain Eqs. (63) and (64) it has been assumed that conditions are such that ‘  L and Eq. (62) is satisfied. Eq. (63) shows that the macroscopic diffusion coefficient is much smaller than the molecular diffusion coefficient D, and can be made arbitrarily small by increasing the potential in the region 0 6 z 6 ‘. By contrast, because the highly charged region of the pore is short compared to the uncharged region, the effect on the partition coefficient is negligible.

[16] [17] [18] [19] [20] [21] [22] [23]

W.M. Deen, AIChE J. 9 (1987) 1409. P. Dechadilok, W.M. Deen, Ind. Eng. Chem. Res. 45 (2006) 6953. P. Dechadilok, W.M. Deen, J. Membrane Sci 336 (2009) 7. G.I. Taylor, Proc. R. Soc. London A 219 (1953) 186. R.M. Weinheimer, D.F. Evans, E.L. Cussler, J. Colloid Interface Sci. 80 (1981) 357. R. Aris, Proc. R. Soc. London A 235 (1956) 67. H. Brenner, Philos. Trans. R. Soc. London 297 (1980) 81. H. Brenner, L.J. Gaydos, J. Colloid Interface Sci. 58 (1977) 312. H. Brenner, P.M. Adler, Philos. Trans. R. Soc. London, A 307 (1982) 149. R.J. Phillips, W.M. Deen, J.F. Brady, AIChE J. 35 (1989) 1761. R.J. Phillips, W.M. Deen, J.F. Brady, J. Colloid Interface Sci. 139 (1990) 363. E.M. Johnson, D.A. Berk, R.K. Jain, W.M. Deen, Biophys. J. 70 (1996) 1017. D.S. Clague, R.J. Phillips, Phys. Fluids 8 (1996) 1720. R.J. Phillips, Biophys. J. 79 (2000) 3350. E.M. Johnson, D.A. Berk, R.K. Jain, W.M. Deen, Biophys. J. 68 (1995) 1561. N. Fatin-Rouge, A. Milon, J. Buffle, R.R. Goulet, A. Tessier, J. Phys. Chem. B 107 (2003) 12126. D. Koch, R.G. Cox, H. Brenner, J.F. Brady, J. Fluid Mech. 200 (1989) 173. W.B. Russel, D.A. Saville, W.R. Schowalter, Colloidal Dispersions, Cambridge University Press, Cambridge, 1989. D.S. Tsai, W. Strieder, Chem. Eng. Commun. 40 (1986) 207. M.M. Tomadakis, S.V. Sotirchos, J. Chem. Phys. 98 (1993) 616. L. Johansson, U. Skantze, J.E. Lofroth, J. Chem. Phys. 98 (1993) 7471. R.J. Phillips, in: J.-P. Hsu (Ed.), Interfacial Forces and Fields, vol. 85, Marcel Dekker, Inc., 1999. R.J. Hunter, Foundations of Colloid Science, vol. I, Oxford University Press, New York, 1989.

260

R.J. Phillips / Journal of Colloid and Interface Science 338 (2009) 250–260

[24] H. Zhou, S.B. Chen, Phys. Rev. E 79 (2009) 021801. [25] P.A. Netz, T. Dorfmuller, J. Chem. Phys. 107 (1997) 9221. [26] T. Miyata, A. Endo, T. Ohmori, M. Nakaiwa, M. Kendo, K.-I. Kurumada, M. Tanigaki, J. Chem. Eng. Jpn. 35 (2002) 640.

[27] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes, Cambridge University Press, New York, 1989. [28] G.D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Oxford University Press, New York, 1985.