8492
J. Phys. Chem. B 1998, 102, 8492-8497
Electrical Interaction of a System Containing Multiple Charged Rigid Spherical Particles Jyh-Ping Hsu* and Bo-Tau Liu Department of Chemical Engineering, National Taiwan UniVersity, Taipei, Taiwan 10617, ROC ReceiVed: January 27, 1998
The electrical interaction of a system containing multiple, charged, rigid spherical particles of arbitrary orientation and size is analyzed theoretically. The equation governing the spatial variation of electrical potential is solved by a boundary integral method, and the electrical free energy of the system under consideration is evaluated. Several approximate analytical expressions for the electrical free energy are presented. We show that pairwise addition will, in general, overestimate the electrical free energy, except that when the electrical double layer around a particle is thin, it is approximately correct. Since the electrolyte concentration at which coagulation occurs is usually high, the classic analysis based on the interaction between two particles is representative of a system containing multiple particles and leads to reasonable results.
Introduction In the analysis of the behavior of a system containing multiple entities it is usually assumed that the interaction between two entities is representative of that of the system. Apparently, this is satisfactory if the concentration of entities is sufficiently low. The spatial variation of the electrical potential for the case of a colloidal suspension, for example, has been analyzed extensively in the literature on the basis of two charged particles.1 In contrast, a direct analysis based on multiple particles is limited. One of the possible approaches for this type of problem is the cell model.2,3 Although it is intuitive and straightforward, this approach is most appropriate for particles having a regular spatial orientation. Hsu and Tseng4 examined the problem of multiple charged entities in an electrolyte solution for the case where the electrical potential is low. An iterative procedure is proposed for the resolution of the linearized Poisson-Boltzmann equation governing the spatial variation of the electrical potential. The criteria of pairwise addition and linear superposition, which are often assumed for a system containing multiple particles, were discussed. While these assumptions may be realistic for point charges and primitive models,5-9 they are not satisfactory for problems related to colloidal particles in which the nature of the charged surface plays a significant role. One of the powerful techniques for the resolution of a boundary-value problem is the boundary integral method,10 in which Green’s theorem is applied to construct an integral solution to a linear differential equation. As the solution is in the form of surface integrals, one need only make numerical evaluations for the domain surfaces and not for the domain interior, which reduces the size of the problem greatly from a three-dimensional to a two-dimensional problem. Or, if there is axial symmetry, it reduces from a two-dimensional problem to a one-dimensional problem. Although Green’s theorem holds only for linear equations, solving the linear form of the equation is an important first step toward solving the full nonlinear form. In addition, it has been found that the linearized form of the Poisson-Boltzmann equation often produces accurate results for electrostatic potentials considerably greater than the linear * To whom correspondence should be addressed. Fax: 886-2-23623040. E-mail:
[email protected].
Figure 1. Schematic representation of the problem under consideration. Rij is the scaled center-to-center distance between particles i and j, ri is the scaled distance measured from the center of particle i, rji is the scaled distance between the center of particle i to a point at the surface of particle j, and aj is the scaled radius of particle j.
assumption limit,11 and from a statistical mechanical basis it may in fact be better justified than the nonlinear form.12 McCarteney and Levine13 and Bell et al.14 successfully applied the boundary integral method in solving the problem of two particles; some approximate analytical results were obtained. The analysis was extended by Sader et al.,15 and a more accurate approximate solution was derived. In the present study, the electrical interaction of a system containing multiple charged particles in an electrolyte solution is evaluated, and a boundary integral method is adopted for the resolution of the governing equations. Analysis Referring to Figure 1, we consider a system containing multiple charged spherical particles in an electrolyte solution. The particles may have different sizes. Let Rij be the scaled distance between the centers of particles i and j, ri the scaled distance measured from the center of particle i, rji the scaled distance from the center of particle i to the surface of particle j, and aj the scaled radius of particle j. These lengths are scaled by the Debye length. Suppose that the electrical potential of the system, ψ, is low. Then its spatial variation at equilibrium can be described approximately by the linearized Poisson-
10.1021/jp9809080 CCC: $15.00 © 1998 American Chemical Society Published on Web 10/03/1998
Electrical Interaction of Multiple Charged Particles
J. Phys. Chem. B, Vol. 102, No. 43, 1998 8493
Boltzmann equation
where
∆y ) y
(1)
κ being the where ∆ denotes the Laplace operator scaled by reciprocal Debye length, and the scaled electrical potential y is defined by y ) eψ/kBT, e, kB, and T being, respectively, the elementary charge, the Boltzmann constant, and the absolute temperature. The reciprocal Debye length κ2 is defined by κ2,
κ2 ) 2e2I/kBT
(6a)
σ ) [σ1,1σ1,2...σ1,Mσ2,1...σN,M]T
(6b)
[ [
χ1,1 ... χ1,N χ) l l χN,1 ... χN,N
(1a)
Gi ) exp(-Fi)/4πFi
(2)
where Fi is the scaled distance from the point charge. Gi is known as the free Green’s function.16 For a system containing N charged particles, the electrical potential can be expressed as N
∫S [Giσi + yiSχi] dSi ∑ i)1 χi ) ∂Gi/∂n
In these expressions, Si denotes the surface of particle i, n the outward normal vector of Si to the bulk phase, and σi the scaled surface charge of particle i defined as
σi ) eσi′/κkBT
(3b)
where σi′ is the surface charge of particle i. We have implicitly assumed that the dielectric constant of a particle is much smaller than that of solution, as is usually the case in practice. It should be pointed out that eq 3 is an exact expression and not based on linear superposition.16 For a charged sphere the singular points of the free Green’s function are located at the particle surface. Therefore, the electrical potential at the surface of particle j can be expressed as
2
yjS )
∫S [Giσi + yiSχi] dSi, ∑ i)1
j ) 1, ..., N
1 l,i ; i ) j and l ) k; i, j ) 1, ..., M; ) - wjχk,j 2 l, k ) 1, ..., N (6f) Gl,k ) [Gi,j,l,k]; i, j ) 1, ..., M
(6g)
l,i Gi,j,l,k ) wjGk,j ; i, j ) 1, ..., M; l, k ) 1, ..., N
(6h)
In these expressions, yi,j and σi,j denote, respectively, the scaled surface potential and scaled surface charge at point j on l,i particle i. χl,i k,j and Gk,j are, respectively, the values of χ and G in which the distance F is measured from the point j of particle l,i k to the point i of particle l. If l ) k and i ) j, then χl,i k,j and Gk,j become singular. In this case, we suggest using the following approach to evaluate their values. For a point at the surface of particle i, we have
∫S Gi dSi ) sinh(ai) exp(-ai) ≡ βii
The value of Gk,j k,j can be evaluated by subtracting the contribution of all other Gl,i k,j from βii that is, M
Gj,j,k,k ) βii -
This is an implicit expression for (or σj), which has the form of an integral equation with an inseparable kernel, and therefore, solving it analytically for yjS (or σj) becomes almost impossible. This difficulty can be circumvented by considering the following approximation for the surface integration of a function F M
(5)
i)1
where wi is a weighting factor and Fi the value of F at point i. The values of wi and M depend on the integration scheme adopted. On the basis of eq 5, eq 4 can be approximated by
χy ) Gσ
k,j wiGk,i ∑ i)1
(8)
i*j
Similarly, we have
yjS
∫F dS ≈ ∑wiFi
(7)
i
(4)
i
(6e)
l, k ) 1, ..., N
N
1
(6d)
(MN)×(MN)
l,i χi,j.l,k ) -wjχk,j ; i * j or l * k; i, j ) 1, ..., M;
(3)
(3a)
(6c)
(NM)×(NM)
χl,k ) [χi,j.l,k]; i, j ) 1, ..., M
i
with
] ]
G1,1 ... G1,N G) l l GN,1 ... GN,N
where I and are the ionic strength and the dielectric constant, respectively. To simplify the mathematical treatment, all the coordinate variables are scaled by κ. For a point charge i, it can be shown by direct substitution that a solution to eq 1 is
y)
y ) [y1,1y1,2...y1,My2,1...yN,M]T
(6)
[ ( )
∫S χi dSi ) - 21 a1i i
1+
]
1 exp(-2ai) ≡ γii ai
(9)
and, therefore, M
χj,j,k,k ) /2 - γii + 1
k,j wiχk,i ∑ i)1
(10)
i*j
Note that G is a symmetric matrix with Gi,j ) Gj,iT, but χ is asymmetric, in general. Equation 6 is a general expression in that an arbitrary surface condition can be assumed. For instance, the surface of a particle can be assumed to have a constant potential, constant charge, or a mixed condition. For illustration, however, we will assume in this discussion that the surface of
8494 J. Phys. Chem. B, Vol. 102, No. 43, 1998 a particle has a constant and uniform potential. The surface potentials of different particles, however, are not necessarily the same. Consider the interaction between a point at the surface of particle j and an arbitrary point at the surface of particle i. Equations 7 and 9 become, respectively,
N
∫S [Gj∆σj] dSj ) ∑ j
i)1 i*j
{
Hsu and Liu
-
βji yiS βii
}
∫S [Gi∆σi] dSi i
(19)
Applying eq 5, we obtain N
ai G dSi ) sinh(ai) exp(-rji) ≡ βji Si i rji
∫
Gj,j∆σj ) -Bj -
(11)
(20)
i*j
and
or
a
[
∫S χi dSi ) 2riji 1 - a1i + i
]
( ) 1+
1 exp(-2ai) ai exp[-(rji - ai)] ≡ γji (12)
∆σj ) -Gj,j-1Bj -
yjS
)
2
{∫S [Giσi] dSi + ∑ i)1 i
yiSγji}
(13)
∆σj ) [∆σj,1...∆σj,M]T
(20b)
Bj ) [Bj,1...Bj,M]T
(20c)
Suppose that σi can be decomposed as
N
σi ) σ0i + ∆σi
N
{∫S [Gi∆σi] dSi + σ0i βji + yiSγji} ∑ i)1
(15)
i
For an isolated particle, the first term on the right-hand side of this expression vanishes. We have
(1/2) - γii i yS σ0i ) βii
N
{
βji
}
∫S [Gi∆σi] dSi + (1/2 - γii)β yiS - yiSγji ∑ i)1 i
ii
(17)
Once ∆σi is determined from this expression, the scaled electrical free energy of the system under consideration, V, with the scaling factor [(κ/)(e/kBT)2], can be evaluated from N
∫S yiS∆σi dSi ∑ i)1
V ) -1/2
(18)
i
Iterative Method Solving eq 17 analytically is almost impossible, and a numerical scheme is preferable, in general. The approach used in the resolution of eq 4 can be adopted here. If the number of particles contained in a system is large, however, the numerical computations involve manipulation of large-scale matrixes. A large memory space is required, and the procedure is time consuming. Also, an ill-defined matrix will create additional problems. To circumvent these difficulties, we suggest using the following iterative procedure, in which eq 17 can be rewritten as
(20d)
ii
In these expressions, βji,k denotes the value of βji evaluated at the point k of particle j. Let ∆σj(m) be the value of ∆σj at the m-th stage. Then ∆σj(m+1) is evaluated by substituting ∆σj(m) into the right-hand side of eq 20a, that is
∆σj(m+1) ) -Gj,j-1Bj -
N
Gj,j-1Gj,i∆σi(m) ∑ i)1
(21)
i*j
We suggest using the initial guess below
∆σj(0) ) -Gj,j-1Bj
(16)
Substituting this expression into eq 15 yields
/2yjS )
βji,k yiS
∑ i)1 β i*j
where is the surface charge for an isolated particle i and ∆σi a perturbed surface charge due to presence of other particles. Substituting eq 14 into eq 13 yields
/2yjS )
Bj,k )
(14)
σ0i
1
(20a)
i*j
N
1
N
Gj,j-1Gj,i∆σi ∑ i)1
where
Since yiS is constant and uniform, eq 4 reduces to
1
Gj,i∆σi ∑ i)1
(21a)
The present iterative procedure is similar to the method of reflections;17 the reasoning, however, is not the same. At each stage of the iteration scheme, the effects of all other particles on the surface charge of a target particle is evaluated. After each particle has been chosen as the target particle once, the surface charge of each particle is adjusted by taking the effects due to the presence of other particles into account. Solving eq 17 iteratively can be avoided under some circumstances. Let us consider two special cases, in which an approximate method can be employed. Method I. Suppose that the separation distance between two particles is large. In this case since the additional charge on particle j, ∆σj, does not affect the additional charge on i, ∆σi, which is due to charges ∆σ0j on all other particles, j, the second term on the right-hand side of eq 19 is negligible. We have
∫
N
[Gj∆σj] dSj = S j
∑ i)1 i*j
{ } -
βji yiS βii
(22)
According to eq 2, Gj is extremely large at a singular point, that is, the integral on the left-hand side of eq 22 can be approximated by integration of the product Gj∆σj near the singular point of Gj. In this case, the ∆σj on the left-hand side of eq 22 can be moved out from the integral sign,13 and we
Electrical Interaction of Multiple Charged Particles
J. Phys. Chem. B, Vol. 102, No. 43, 1998 8495
have N
∆σj ) -
βji
yiS ∑ i)1 β β i*j
(23)
ii jj
Substituting this expression into eq 18 gives N-1
V ) 4π
N
aiajyiS yjS ∑ ∑ j)1 i)j+1
exp[-(Rij - ai - aj)] Rij
(24)
This is equivalent to saying that the total electrical free energy can be obtained by summing all the paired free energies between two particles. Here, the paired free energy has the Yukawa form, which is derived on the assumption that the electrical potential can be obtained through linear superposition. Method II. Suppose that the double layer around a particle is thin. By applying the same argument as that used in the previous method, ∆σi in the first term on the right-hand side of eq 17 is moved out from the integral sign, and we have N
N
βji
∆σiβji ) - ∑ yiS ≡ Bj; ∑ i)1 i)1 β i*j
j ) 1, ..., N
(25)
ii
The solution to this equation can be expressed as
[ ] [
∆σ1 β11 ... β1N l )- l l ∆σN βN1 ... βNN
][ ] -1
B1 l BN
(26)
It should be pointed out that each component of the matrix on the left-hand side of this expression is a function. For example, if N ) 2, then eq 26 leads to
∆σ1 ) -
β22B1 - β12B2 β11β22 - β12β21
(27a)
∆σ2 ) -
β11B2 - β21B1 β11β22 - β12β21
(27b)
In the case where particles are regularly arranged, eq 26 leads to N
∆σ1 ) -B1/
β1i ∑ i)1
(28)
Discussion For illustration, we assume that the system is monodispersed in the numerical simulation. The weighting factor wi in eq 5 is chosen on the basis of the Simpson’s 1/3 rule. Figure 2 shows the variation of the scaled electrical free energy, V/(Ny2S), as a function of the ratio (scaled center-to-center distance between two particles/scaled particle radius), (R/a), for the case of two particles at various a, and Figures 3 and 4 show the results for the cases of three and four particles. For a fixed particle size, the greater the value of a, the thinner the electrical double layer around a particle. Both the results based on the present boundary integral method and those based on the approximate methods (methods I and II) are shown in Figures 2-4; the results based on pairwise addition are illustrated in Figures 3 and 4 for comparison. Figures 2-4 reveal that the greater the number of particles contained in a system, the greater the electrical free
Figure 2. Variation of the scaled electrical free energy, V/(Ny2S), as a function of the ratio (scaled center-to-center distance between two particles/scaled particle radius), (R/a), for the case N ) 2. Key: s, present boundary integral method; O, approximate method I; ×, approximate method II. Parameter used: (a) a ) 0.1; (b) a ) 1; (c) a ) 5. All particles are identical with the same scaled radius a.
8496 J. Phys. Chem. B, Vol. 102, No. 43, 1998
Hsu and Liu
Figure 3. Variation of the scaled electrical free energy, V/(Ny2S), as a function of ratio (scaled center-to-center distance between two particles/ scaled particle radius), (R/a), for the case N ) 3. Key: same as Figure 2 plus ∆, pairwise addition.
Figure 4. Variation of the scaled electrical free energy, V/(Ny2S), as a function of the ratio (scaled center-to-center distance between two particles/scaled particle radius), (R/a), for the case N ) 4. Key: same as Figure 2 plus ∆, result based on pairwise addition.
energy, and the greater the (R/a), the smaller the deviation caused by using the methods I and II and pairwise addition. These figures also reveal that the greater the a, the smaller the
deviation caused by using the approximate methods and pariwise addition. This implies that the thinner the electrical double layer, the better the performance of these approximate methods.
Electrical Interaction of Multiple Charged Particles Figures 3 and 4 suggest that the performance of these methods, from the worst to the best, follows the order (method I) < (pairwise addition) < (method II). That of method II is satisfactory, except when (R/a) is small, that is, when the particles are close together. As revealed by Figures 3 and 4, pairwise addition overestimates the electrical free energy, in general, except when the electrical double layer around a particle is thin. This is consistent with the result obtained for the case of regularly arranged infinite cylinders.2 For the case where the electrical double layer is thin, pairwise addition is appropriate, since each particle behaves like an isolated entity. In practice, the electrolyte concentration at which coagulation between particles occurs is usually appreciable; that is, the electrical double layer is thin. This is why the classic analysis based on the interaction between two particles is representative of a system containing multiple particles and leads to reasonable results. On the other hand, if the electrical double layer is thick (or a is small) and the particle surface is maintained at constant potential, the resistance experienced by one of three particles as it approaches the other two particles is about the same order as when it approaches a single, but bigger, particle having the same surface potential. In this case the electrical free energy for three particles is about the same order as that between a particle and a bigger one. The latter, however, is smaller than the electrical free energy of two pairs of smaller particles. This is why pairwise addition overestimates the electrical free energy. As compared to the boundary collocation method, determining the form of the eigenfunction of eq 1 is unnecessary in the present boundary integral approach. This is nontrivial in the former, especially when multiple particles are present, since defining a global coordinate system becomes a difficult task.
J. Phys. Chem. B, Vol. 102, No. 43, 1998 8497 Another disadvantage of the boundary collocation method is that the boundary conditions associated with eq 1 are satisfied only at the collocation points specified, and appreciable deviations may be observed at other points. In contrast, this does not occur in the present boundary integral method. This is because the discretized expressions, such as eqs 6 and 20, which are based on eq 5, are sufficiently accurate provided that the functions representing the boundary conditions are smooth. References and Notes (1) Hunter, R. J. Foundations of Colloid Science; Oxford University: London, 1989; Vol. 1. (2) Reiner, E. S.; Radke, C. J. Am. Inst. Chem. Eng. J. 1991, 37, 805. (3) Van Keulen, J.; Smith, J. A. M. J. Colloid Interface Sci. 1995, 170, 134. (4) Hsu, J. P.; Tseng, M. T. J. Colloid Interface Sci. 1996, 184, 289. (5) Chu, X.; Wasan, D. T. J. Colloid Interface Sci. 1996, 184, 268. (6) Blum, L. Mol. Phys. 1975, 30, 1529. (7) Blum, L.; Hoye, J. S. J. Phys. Chem. 1977, 81, 1311. (8) Waisman, E.; Lebowitz, J. L. J. Chem. Phys. 1972, 56, 3086. (9) Gillan, M.; Larsen, B.; Tosi, M. P.; March, N. H. J. Phys. C: Solid State Phys. 1976, 9, 889. (10) Brebbia, C. A. The Boundary Element Method for Enginners; Willey: New York, 1978. (11) Fushiki, M.; Svensson, B.; Jonsson, B.; Woodward, C. E. Biopolymers 1991, 31, 1149. (12) McQuarrie, D. A. Statistical Mechanics; Harper & Row: New York, 1976. (13) McCartney, L. N.; Levine, S. J. Colloid Interface Sci. 1969, 30, 345. (14) Bell, G. M.; Levine, S.; McCartney, L. N. J. Colloid Interface Sci. 1970, 33, 335. (15) Sader, J. E.; Carnie, S. L.; Chan, D. Y. C. J. Colloid Interface Sci. 1995, 171, 46. (16) Greenberg, M. D. Application of Green’s Functions in Science and Engineering; Prentice Hall: Englewood Cliffs, NJ, 1971. (17) Happel, J.; Brenner, H. Low Reynolds Number Hydrodynamics: with Special Applications to Particulate Media; Nijhoff: Boston, 1983.
