- Email: [email protected]
Numerical solution of the Poisson—Boltzmann equation
Numerical Solution of the Poisson-Boltzmann Equation A. E. JAMES AND D. J. A. WILLIAMS* Department of Oceanography and *Department of Chemical Engineering, University College of Swansea, Singleton Park, Swansea, United Kingdom
Received March 8, 1984; acceptedJanuary 2, 1985 A Galerkin finite element scheme is used to provide flexible numerical solutions to the PoissonBoltzmann equation for electrical double layers in one and two dimensions. A Newton sequence is used for the solution of the set of nonlinear equations arising from the finite element discretization procedure. Convergence of the Newton sequence is rapid, generally occurring in less than four iterations. Isoparametricfinite elements are used for computational ease. Error analysis is applied to the solution of the linearized Poisson-Boltzmann equation (Debye-Hiickelapproximation) and also to the full iterative scheme. Numerical solutions are presented for (i) overlapping cylindrical and planar electrical double layers, and (ii) the overlapping electrical double layers of two interacting platelets. © 1985 Academic Press, Inc. INTRODUCTION
colloidal particles. Analytic solutions to the Poisson-Boltzmann equation are limited to special cases in one dimension. Complete solutions for the electrical potential decay from an infinite flat plate and have been given by Gouy (4) and Chapman (5) who consider the plate to be immersed in a binary symmetrical electrolyte. Nir and Bentz (6) and Bentz (7) have considered the effects of mixtures of electrolytes. The potential distribution for other geometries has been given for single spheres (Duchin et al. (8); Sigal and Shamanski (9); Stigter (10); Brenner and Roberts (11); Abraham-Schrauner (12) and White (13)) and single cylinders (MacGillivray and Winklemann (14); Duchin et al. (8); Sigal and Shamanski (9); Mills (15); Philip and Wooding (16); Stigter (17); van der Drift (18); Torrie and Valleau (19)). When the surface potential is 425 mV it is possible to linearize the exponential term in the Poisson-Boltzmann equation. This simplification is known as the Debye-Hiickel linearization and solutions to this equation are termed Debye-Hfickel solutions. The linearization removes the dependence of the solution on specific electrolyte types, so that Debye-Hiickel solutions are readily applica-
A knowledge of the distribution of electrical potential about colloidal particles dispersed in electrolyte is of great practical and theoretical importance. Typically, it is needed to calculate the free energy of interacting electrical double layers, which together with the attractive van der Waal's energy and energy attributable to hydrodynamic interactions, determine the stability of colloidal systems. The notion of the electrical double layer was proposed by Helmholtz (1), who likened it to a simple electrical capacitor. Electrical double layer theory has been modified by Stern (2) and Grahame (3). These authors split the electrical double layer into two regions: the inner or fixed region and the outer or diffuse region. The former is determined by the closest approach of counterions and the variation of potential in the region is likened to that in an electrical capacitor, In the diffuse layer, the distribution of potential is described by the Poisson-Boltzmann equation. The Poisson-Boltzmann equation is the fundamental differential equation describing the distribution of electrical potential about 44 0021-9797/85 $3.00 Copyright © 1985 by Academic Press, Inc. All rights of reproduction in any form reserved.
Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION
ble to electrolyte mixtures. Debye and Hiickel (20) give a solution for a single spherical electrical double layer, while Dube (21) and Sparnaay (22) give results for a single cylindrical double layer and two overlapping cylindrical double layers, respectively. The case of two parallel plates has been analysed by Hogg, Healy, and Fuerstenau (23). Numerical solutions to the Poisson-Boltzmann equation have been given for a single sphere (Muller (24); Guggenheim (25), and Loeb, Overbeek and Wiersema (26)), two identical spheres in a binary symmetrical electrolyte (Hoskin (27); McCartney and Levine (28)), a single cylinder (Stigter (17)), and the potential distribution between two parallel plates in mixed ionic systems (Bresler (29)). Analytic solutions to the Poisson-Boltzmann equation are limited to regularly shaped individual colloidal particles. Similar restrictions also apply to existing numerical solutions. Numerical solutions of the PoissonBoltzmann equation have used the finite difference method combined with first-order iterative methods and have not been flexible enough to allow for changes in size, shape, and geometry of the system under consideration. Currently most analytic and numerical solutions to the Poisson-Boltzmann equation have been one-dimensional simplifications which have necessarily limited application in real situations. There is need for a general numerical solution to the Poisson-Boltzmann equation for arbitrary geometry. To this end, a finite element scheme is presented using a second-order iterative scheme to speed the convergence of the nonlinear set of equations obtained from the discretization process. 2. THEORY
The Poisson-Boltzmann equation written for a two-dimensional electrical double layer is 020
020
~ 2 + 0 Y~ = -0(0) OX
[ 1]
where X and Y are nondimensional distances,
45
0 is the nondimensional electrical potential, and p(O) is the nondimensional space charge density. That is X = rx, Y = ry, and ¢ = op/ kT, where x and y are Cartesian coordinate distances, K = (2NAIe2/~kT) 1/2, e = electronic charge, ~b = electrical double layer potential, k = the Boltzmann constant, T = absolute temperature, ~ = permittivity of the electrolyte, NA = Avogadro's number, I = (ZV2hCh)/2 = the ionic strength of the electrolyte, with vh and ch being the valency and concentration of the hth ionic species. The summation is carried out over all ionic species in the electrolyte and vh contains the sign of the charge of the hth ionic species. Consistent S.I. units have been used. For example the concentrations ch and I are in kmole m -3 and K the DebyeHfickel parameter has units of m -1. The nondimensional space charge density is written
p(O) = (r--g---ekT) ~ VheNAche-"h* 1
- - - Z VhChe-vh*. 21
[2]
The summation again is taken over all ionic species present in the electrolyte. In order to solve Eq. [ 1] it is necessary to define the potential, 0, and the potential gradient 00/0N at the boundaries, 1", of the chosen solution domain, ft. Here N (=Kn) is the outward normal direction from the boundary 1". The boundary condition is written as ¢ = 0b on r I [3] and 00
--KO"
ON - 2NAIe on 1"2
[4]
where F1 includes all parts of F where the nondimensional potential is fixed, 1,2 includes all parts of F where the nondimensional potential gradient is defined, and a is the surface charge density (C m -E) on 1"2, (ga/ 2NAIe) is the nondimensional surface charge density. In order to produce a numerical solution to the Poisson-Boltzmann equation using Journal of Colloid and Interface Science,
Vol. 107,No. 1, September1985
46
JAMES A N D W I L L I A M S
the finite element technique it is necessary to divide the solution domain into a grid (as in the finite difference procedure). This grid need not be regular and the finite elements thus produced are usually triangular or rectangular. In the present work the latter type were used. Each finite element has a fixed number of nodal points along its boundary and the numerical solution is generated at these points. The variation of a parameter (e.g., 4, X, Y) within a finite element is described in terms of the nodal values. Thus for a given area, the accuracy of the approximation increases as the number of finite elements increases and their size decreases. The finite element method is an extension of the Rayleigh-Ritz-Galerkin technique, whereby the functional form of the original differential equation is discretized by assuming the solution to be given by sets of overlapping polynomials whose nodal values have to be determined. The discretization process is applied to each finite element in the solution domain so that the variation of the unknown variable in each finite element is described by a set of polynomials which are known as shape functions. The order of the polynomial used characterizes the type of finite element by specifying the number of nodal points on the side of a finite element. The discretized set of equations describing the variation of the unknown parameters in individual finite elements is assembled using the direct stiffness method (Desai and Abel (30)) into an overall or global set, which is subsequently solved for the nodal values of the unknown variables. In the present work the Galerkin (31) method is used to generate the equivalent functional form of the Poisson-Boltzmann equation, and this is subsequently discretized. The usefulness of this method was recognized by Oden (32, 33) and Szabo and Lee (34); further, Oden (35) suggests that it is the most powerful technique for generating finite element representations of nonlinear differential equations. The Galerkin method does not operate directly with the relevant differential Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
equation but with the equivalent functional form of the continuous problem, of which the solution is assumed to be a combination of trial functions, say Pi- Thus at any node i within a finite element having a total of m nodes it follows that the equivalent functional form of Eqs. [1], [3], and [4] is y,/d2q5
1
02~b
VhChe_Vh+)
f~e
× d~2e + f~ Pi(X, Y)(c~ -
4)b)dF~
dr
I Ka
04,]dF~
wherein Pi(X, Y) is the weighting function with respect to the ith node and the superscript e refers to a single finite element. Green's lemma (see Wylie (36)) is used to convert the second-order partial derivatives in Eq. [5] into products of first-order terms (see Appendix II for details). The transformed equation is written
~e
,ax o x
& 2I (~ VhChe-v*~)dXdg Kff
e
+ ~r Pi(2-~Aie)dr2=O.
[6]
It is worthwhile noting that the second term in Eq. [5] vanishes as (4) - 4)b = 0) by definition on r~, while the third term in Eq. [6] is the line integral of the nondimensional surface charge along F~. The first two terms are evaluated over the area of the finite element and implicitly contain boundary conditions of the form 4~ = q~b since this condition is specified at any node on the boundary. The final step in the finite element formulation is the discretization of the equivalent functional form of the Poisson-Boltzmann equation given by Eq. [6]. Within any
47
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION
particular finite element the variables X, Y, and ¢ can be written in terms of their nodal values as x e : E pe(~, n)X e = [e]e{x}e
ye = X pe(~, y)y~
=
[p]e{y}e
q~e = E pe(~, n)¢e = [p]e{¢}e.
[7]
charge, % along each side of a rectangular finite element in turn. This is consistent with the convention that the solution domain is always to the left of the line along which the integral is performed. The values A, (i = 1, 2, 3, 4), are given by Eqs. [B17]-[B19] while the Jacobian determinant IJI --- O(S, Y)/O(~, n). These parameters are required to change from the cartesian (X, Y) coordinate system to the curvilinear (~, n) coordinate system in which the integration of the equivalent functional form of the Poisson-Boltzmann equation is performed. Equation [8] can conveniently be incorporated into a matrix equation for the m nodes in a finite element [F]e{dp}e+ {C(~b)}e-[ - { S } e = { 0 } e [9]
The summations are made over the m nodes in the finite elements; [p]e is the (1 × m) line vector of shape functions; {X}, {Y}, and {¢} are the (m × 1) column vectors of the nodal values of X, ]1, and ¢, respectively. P~((, n) is a polynomial of chosen degree (usually less than 3) and ( and n are curvilinear coordinates introduced to enable the integration to be performed over a uniform region. If the same polynomials, Pg, are used where [F]e is an (m × m) symmetric positive for both the weighting function in Eq. [5] definite matrix of which the general member and also for the description of a variable f0 is written within a finite element then the finite elements are isoparametric (Zienkiewicz (37)). f/j = :__'-11__'-I 1 ~ 1 Equation [7] is substituted into Eq. [6] to give (OPiOPj OPiOPJ\jd d
×\OXOX+-~-~)l
f~'-I f~"l (Opi O[p]e Opio[p]e ) 1 l ~OX OX {l#}e"~ OY OY {~)}e
[101
The (m × 1) column vector {C(¢)} contains the nonlinear terms related to the space charge density and has the form at the ith row of
f~'-lf-lpi 1 I 2I
× lJld~dn -
[ ~ 7/.
× (X VhChexp(-- Vh[Ple{¢}e))lJldtidrl cR¢)
= 1
-[-
Pi E (eql~ ql)mld~ 1 ql
+
£1Pi(~ Pq2~q2)A2dn
+
j
1
1
× ( Z VhChexp(--Vh[p]e{¢}e))lJld~dvl [1 1] and the surface charge boundary condition vector, (m × 1) has the general term
q2
~- I +1
L
Se=
Pi(Z Pq3~q3)A3d~
Pi(~_,Pql~ql)Ald~
I
q3
ql
j
"+ 1
+ +
Pi(E Pe---e)A4dn +1 q4
= 0
[81
where ~q = (K/2NAIe)~q, and ---q = 0 if q is not a boundary condition node. The summations ~q~, ~]q2, ~q3, and Zq4 refer to the evaluation of the nondimensional surface
-1
+ J+l
Pi( ~ Pq2==-q2)~2d~ q2
Pi(~ eq3~q3)m3d~ q3
+
f? P,(Z Pe%4),~dn.
[12]
1 q4 Journalof ColloidandInterfaceScience,Vol. 107, No. 1, September1985
48
JAMES AND WILLIAMS
This last term assumes nonzero values only at nodes where surface charge is specified as a boundary condition. Because of the nonlinearity of terms in the space charge vector Eq. [9] cannot be inverted immediately. In the present work, a Newton sequence (see Rail (45)) is used to solve the nonlinear set of equations described by Eq. [9]. Let a nonlinear operator A~(~) be defined as
A~(4~) = [Fle{~b} + {C(q~)} eq-
{S} e
[131
It is apparent from Eq. [18] that [ J A G ] e is symmetrical, which is useful in the numerical computations. The actual Newton sequence required is now written as
[SAC]e{ Adp}e*+l = - [ F y { q ~ } e . - { C ( 4 , ) } e . - {S}~n. {4~}~n.+1 = {~}e. + {dX4~}e.+l
[19]
3. COMPUTATIONAL ASPECTS OF THE SOLUTION
SO that the Jacobian matrix A~(q~) is written The Poisson-Boltzmann equation has been discretized and a suitable Newton sequence OA~ck Ael(dP) - - - [ J A C ] e. [141 has been presented for the solution of the resultant set of nonlinear equations describing the potential distribution within a single finite The Newton sequence required to solve element. The solution for the potential A~(qS) = 0 is defined by throughout an entire region which has been A~(~b){A~b}e*+, = -A~(4~),* [15] divided into finite elements is obtained by constructing global matrices corresponding {q~}e.+l = {q~}ne. q- {Aq~e}n,+l to [JAC] e and [F] e and the global column where the subscript n* refers to values of ~b vectors {C(q~)}e and { S } e. This global assemand 2x~ at the n*th iteration. An initial set bly is achieved using the direct stiffness of values to start the sequence is taken from method (e.g., see Oden (35) or Desai and the solution to the Debye-Hiickel problem Abel (30)). A frontal solution scheme (see (see Appendix I). The convergence of Newton Irons (38)) is used to invert the global set of sequences is discussed in detail by Rail (38). equations. The Newton sequence does not Practically, A~(q~) is an m x m matrix, the require the recalculation of [F] e on each iteration. The coefficients in every other mageneral member being given by trix and vector have to be evaluated at each O(a~(ck)) jac~ - - [ 16] step. Provided a good initial guess is made 0~j using the Debye-Hfickel solution (see ApHere a~i(ck) is the i th row of A~(q~) such that, pendix I) the Newton sequence usually converges within three iterations. Programs imusing Eq. [ 11 ], plementing this numerical scheme have been m written in F O R T R A N for an ICL 1904S aoi = ~ f,~q~j + ce(qs) + s e. [17] computer. as
j=l
Thus, using Eqs. [10], [11], [16], and [t7] m and noting that q5 = ~i=l PigJi, 1 f+l
jac~=f°
e__
~d-i
i'+l
,-i
pipj
p*
x ( Z VhChexp(-vhck))lJId~dn.
[18]
h=l Journal of Colloid and Interface Science, Vol. 107,No. 1, September1985
4. ERROR ANALYSIS
There are no analytic solutions to the Poisson-Boltzmann equation in two dimensions even for such a simple geometry as a single isolated charged cylinder in a binary symmetrical electrolyte. Although this is a one-dimensional problem in cylindrical coordinates, it would have provided a suitable
NUMERICAL
SOLUTION
OF THE POISSON-BOLTZMANN
test in rectangular Cartesian coordinates. There are a number of approximate analytic solutions but these are of little use in ascertaining the accuracy of the finite element solution. In one dimension there is an analytic solution to the Poisson-Boltzmann equation (Gouy (4) and Chapman (5)) for binary symmetrical electrolytes. The equation describing the decay of potential from an infinite fiat plate into an infinite fiat double layer is written as e~°/2 + 1 + (e ~°/2- l)e - x 4 ' = l n e ~ o / 2 + 1 - ( e *°/2- l)e - x "
[20]
Here 4)0 is the reduced surface potential. For 4'0 < 1 the corresponding Debye-Hfickel solution is [21 ]
4) = 4'oe-x.
Equations [20] and [21] were derived using the boundary conditions X=0,
4'=
4,0
X= ~,
~-~=0.
a4'
[22]
The infinite boundary condition is problematical in some numerical solutions. In the present work, it was sufficient to choose a value of X large enough that d4'/dX _~ O. In practice, this occurs at around 30 units from any surface. The boundary conditions used in the flat plate solutions are X=O
4'=4'0
X = 30,
d4'/dX = 0.
49
n is the order of the functional that has been constructed from the original differential equation, and B is a constant. In the present work n = 1, and k = 2, 3, 4, corresponding to linear, quadratic, and cubic finite elements. Equation [24] suggests that the error in each of these finite elements is proportional to h 2, h 4, and h 6, respectively. The errors between the finite element solution to the distribution of electrical potential in a plane double layer and Eqs. [20] and [21] have been computed and are shown in Figs. 1 and 2. The values of B, obtained by regression of ]error[ o n h (2K-n) are shown in Table I. Although the finite element solution to the Poisson-Boltzmann equation uses an iterative scheme the values of B are not significantly different from those seen for the linear Debye-Hfickel solution. It is notable that the Newton sequence takes no more than three iterations to converge while other iterative procedures used to solve the PoissonBoltzmann equation have used over 20 iterations (see Hoskin and Levine (27); McCartney and Levine (28)). 5. E R R O R S
IN TWO-DIMENSIONAL SOLUTIONS
There are no suitable analytic solutions to test the convergence of the numerical solution 0-
-2-
[231
The forcing of the potential gradient at large X does not introduce any significant errors into the computer solution especially as an iterative technique is used. Strang and Fix (39) have shown that the simplest and most fundamental error estimate to linear finite element problems is given by lerrorl --< B2h 2(k-n)
EQUATION
[241
where h is the length of the element side, k is the order of the shape functions (P;((, n)),
L
o
_/,-
-6
-8
-10 -2 log
(h)
FIG. 1. The relationshipbetween [error[and h, length of a finiteelement,when linear(LFES),quadratic (QFES), and cubic(CFES)finiteelementsare used in the numerical solution of the Poisson-BoltzmannEquation. Journal of Colloid and Inte(face Science, Vol. 107, No. 1, September 1985
50
JAMES AND WILLIAMS 0
Figure (3) shows the arrangement of finite elements and relevant boundary conditions -2 for the numerical solution to the PoissonBoltzmann equation about a single cylinder. Typical results for the Debye-HiJckel solution (Eq. [24]) and the numerical solution using -6 quadratic finite elements, are shown in Table o= II. It can be seen in Fig. 3 that for the aspect -8 ratios of the elements to be kept in the region of unity, it is necessary to increase the density of elements (i.e., reduce their size) as the -10 -2 cylinder wall is approached. In practice, it log (h) was found that it was not necessary for the FIG. 2. The relationship between [errort and h, the aspect ratio of elements close to the cylinder length of a finite element, when linear (LFES), quadratic to be close to unity for good numerical (QFES), and cubic (CFES) finite elements are used in the numerical solution of the Poisson-Boltzmann Equa- solutions to be obtained. The main criteria for the simple single cylinder geometry seems tion. to be the radial length of the element, which once small enough, ensures good agreement to the Poisson-Boltzmann equation with re- between analytic (Eq. [25]) and numerical spect to element size for a two-dimensional schemes. This result is possibly linked to the problem. It would seem reasonable that the constraint that for a single cylinder &b/ON relationship between the length of an element = 0, where N is the direction normal to the side and the absolute error will hold in both radial direction. When the geometry becomes Cartesian (X, Y) directions if finite elements more complex, e.g., two or more cylinders, with an aspect ratio of unity are used (ratio Odp/ON d= O. of m a x i m u m length to m i n i m u m length, see Desai and Abel (30)). A major attraction of 6. TWO-DIMENSIONAL NUMERICAL the finite element method is its flexibility in SOLUTIONS allowing elements of different size and shape One-dimensional solutions to the Poissonto be used in the same mesh; clearly, this attraction is lost if all the elements are equi- Boltzmann equation have limited applications sized. Further, in m a n y two-dimensional in colloid chemistry. In m a n y instances it is problems, it is not possible to use rectangular necessary to know the electrical potential elements with an aspect ratio of unity. It can be shown (Sparnaay (22)) that under the restrictive conditions of the D e b y e TABLE I H/ickel approximation, an analytic solution Values of B, Determined by Linear Regression of the can be obtained. That is Error Relationship (Eq. [25] and Fig. 1) for Linear, Qua4~ = 4~oKo(R)/Ko(Ro) [25] dratic, and Cubic Finite Elements
j
where R is the nondimensional radial distance from the axis of the cylinder and R0 is the cylinder radius. ~bo is the reduced surface potential of the cylinder. K0( ) is a modified Bessel function of the second kind and zero order. Journal of Colloid and Interface Science, Vol. 107, No. I, September 1985
Linear element Quadratic element Cubic elements
Debye-Hfickel solution
PoissonBoltzmann solution
5.46 × 10-2 7.03 × 10-3 10-3
7.16 X 10 -2 7.41 X 10-3 3.16 × 10-4
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION
element
51
solution
FIG. 3. Arrangement of finite elements and relevant boundary conditions for a single cylinder.
distribution in the double layers about two or m o r e particles immersed in electrolyte. Some interesting and relevant geometries are those involving cylindrical or plate-like particles. A two-dimensional solution to the P o i s s o n - B o l t z m a n n equation can be utilized in such cases and others where s y m m e t r y can be used to reduce the dimensionality o f the problem. The problem o f interacting particles having different geometries m a y either be a real case where the geometry is known, e.g., between pairs o f tobacco mosaic virus (cylinder-cylinder) (see Brenner and McQuarrie (40); James and Williams (41)), between protein molecules and a cell wall (cylinder-plate), or in a less well defined system, where some geometric shape m u s t be assumed to evaluate interparticle electrical double-layer potentials,
e.g., kaolinite suspensions (Flegmann, G o o d win, and Ottewill (42); Williams and Williams (43); James and Williams (44)). Firstly, the problem o f overlapping cylindrical and planar electrical double layers typifying the protein cell wall interaction is examined. Figure 4 shows the b o u n d a r y condition o f the cylinder plate geometry. The general layout o f the finite element mesh is shown in Fig. 5. A n u m b e r o f mesh refinements were required to obtain a convergent solution. The areas o f modification are shown in Fig. 5. Figures 6 and 7 show results for the electrical double-layer interaction between a cylinder and a plate at different separations. The increase in potential gradient reflecting an increase in the attractive free energy is apparent at the smaller separation. The distortion o f the potential distribution about the
TABLE II Comparison between Debye-Hiickel Solution for a Cylinder and the Analytic Debye-Hiickel Solution (Eq. [25]) for the Potential about a Cylinder, Using Quadratic Finite Elements~ Number of elements R
26
39
49
56
66
Anflytic equation [25]
0.2 0.7 1.2 2.2 4.2 6.2
2.0 0.806 0.380 0.107 1.87 X 10 -2 2.07 X 10 -3
2.0 0.759 0.367 0.103 1.02 X 10 -2 1.15 X 10-3
2.0 0.755 0.364 0.102 1.02 X 10-2 1.15 X 10-3
2.0 0.7543 0.3637 0.1019 1.019 X 10 -2 1.145 X 10-3
2.0 0.7538 0.3634 0.1019 1.019 X 10 -2 1.141X 10-3
2.0 0.7538 0.3634 0.1019 1.019 X 10 2 1.141X 10 -3
a
Cylinder radius = R0 = 0.2. Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
52
JAMES
AND
9¢ - - =0 9N
to describe the interaction and geometry of two interacting platelets, which could be clay particles. Typical solutions are illustrated in Figs. 10 and 11. The potential gradient is seen to increase at the closer separation (Fig. 11). It is also clear that the electrical field associated with the edge becomes more localized and distorted at close separation (cf. Fig. 9). Both of these observations are consistent with an increase in the attractive free energy of the system.
y=÷-~
¢
/
9~ =0 9N X= - e e
:
WILLIAMS
~02
¢=~0
7. S U M M A R Y
9 ~ = 0 along AA' 9N /
S /
9~ =0 9N
--
V = -~ FIG. 4. B o u n d a r y
conditions for the cylinder-plate
geometry.
cylinder shows the effects of an increasingly strong interaction at close separation. Figures 8 and 9 show the boundary conditions and basic finite element mesh used
;
This work has shown that in one or two dimensions, the finite element method is ideally suited to solutions of the nonlinear Poisson-Boltzmann equation. The variable geometry of interacting charged surfaces is more readily accommodated by a flexible, nonuniform finite element mesh than by a finite difference scheme. The finite element procedure produces a set of nonlinear equations which are solved in the present work with the aid of a Newton sequence; this sequence converges rapidly when a trial so-
~N =0
3
3
3
3
3
3
o2
1
First
modification 2
Second modification
3
Third modificotion
+o
9N =0
Ol
--=0 9N
FIG. 5. Finite e l e m e n t m e s h for cylinder-plate geometry.
Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
53
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION
-5.0 -L.O
-30
-2.0
-1.0
0.0
1.0
2.0
3.0
4.0
5,0
6.0
7.0
8.0
9.0
X
FIG. 6. Distribution of potential ~b, for the interaction of a cylinder (Ro = 1.0) with a plate at a separation X = 8.0.
lution obtained from the linearized DebyeHiickel equation is used. Strang and Fix (39) have given a criterion for the convergence of finite elements in linear problems (e.g., the Debye-Hiickel equation). Interestingly the present work has shown that the error relationship is applicable to the nonlinear Poisson-Boltzmann equation, when either linear, quadratic, or cubic finite elements are used. The finite element method is a very convenient and powerful method of studying the various modes of potential interaction that can occur between particles of irregular geometry. In contrast to finite difference schemes the actual geometrical surface of the particles can be easily incorporated in the numerical scheme. The technique is readily extended to three-dimensional problems such
as the potential interactions between pairs of spheres and to more complicated systems, such as arrays of particles. APPENDIX I: FINITE ELEMENT FORMULATION FOR THE DEBYE-HOCKEL APPROXIMATION
The Debye-Hiickel approximation to the two-dimensional Poisson-Boltzmann equation is written 0%
0% + ~ = 4,
[A1]
with boundary conditions ~b = 4~b on P L
Odp ON
[A21
K~
- -
2NAIe
o n 172.
[A3]
Journal of Colloid and Interface Science, Vol. 107, No. 1, Septembert985
54
JAMES A N D WILLIAMS
7.0 6.0 5.0 &0
Note that this equation differs from Eq. [6] in that the second term contains no nonlinear function of 4. Substituting for 4~ and using isoparametric finite elements leads to the matrix equation for a single finite element
3.0
[Dle{dp}e + {S} e =
0,
[A6]
20 1.0 0.0 -19 -2.0
where
+l (Op~opj OPiOPj ) + -if-f-if-f+ PiP/ 4J= f-l f-, tOX OX x IJld~n [A7] +1
-3.0 -/,.0 -5.0
and
Si =
£.-1Pi(Z Pql-~ql)Ald~ +F' 1
-6,0 -70
ql
1 -5.0
-40
-3.0 -2.0
-1.0
00
1,0
20
30
4.0
q2
5.0
X
Pi(
FIG. 7. Distribution of potential 4, for the interaction of a cylinder (R0 = 1.0) with a plate at a separation X = 4.0.
1
Z
Pq3m-q3)A3d~
q3
j
,t_ 1
+
Pi(~ Pq4=--q4)A4dT1• [A8] +1
q4
Galerkin's method is applied to Eq. [A1][A3], so that at any node i in a finite element, it follows that / / / / /
Y =~
f f Pi,ox2 0y2-
9N
¢=~
X=0
+ fr~ P~(4a- 4,b)dr[
~N
ra're 2 = 0. [A4]
+fr~Pi
[oei 04
% 0/
¢~-,~
A'
/
+ ~ ) d
After applying Green's lemma and some rearrangement (see Appendix II), Eq. [A4] is written
ff
A
/ / / /
ori~y)dfU+ ff e''d"e
£e
£e
+
dr2 = 0.
[AS]
Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
/
¢ / / / / /
X = 0 9~
:0
~N
¢ 9~=0 ~N
FIG. 8. Boundary conditions for the interaction of two plates in edge-face mode.
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION
9~
9N
55
-0 02
02
01
~ "~N
0
FIG. 9. Finite element mesh for two interacting plates in edge-face mode.
4.0
5.0
6.0
7.0
8.0
9.0 10.0 11.0 12~ 13.0 l&0 15.0 16.0 1Z0
X FIG. 10. D i s t r i b u t i o n o f potential th, f o r the edge-face interaction o f t w o plates at a separation X = 8.0.
Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
56
JAMES
0.02 0.05 0:1 02
AND
WILLIAMS
Proof of Eq. [B1] may be found in most mathematical texts concerning line, surface, and volume integrals. (For example, see Wylie (36) p. 568 et seq.) On substitution of
05 1C ZO 50 5.0 4.0 3.0
v =e/0-0q~
20 1.0
T
0.0 -1.0
Y
[B21
and
04
[B3]
u = - P i O---Y Eq. [B 1] is written
-20 -3.0 -4.0
£
-5.0
o
o2
r
= fr ( - e i - ~ dX + Pz ox dY) d ,
-6.0
[B41
-ZO
so that 3.0
&0
50
6.0 X
7.0
&0
02, oei0,+
g.0 10.0 11.0 12.0 13.0
=,
F I G . 11. Distribution of potential ~, for the edge-face interaction of two plates at a separation X = 4.0.
£
04
04
× dXdY = fr Pi(ff-x dY - -~ dX ) • The matrix [D] e is symmetrical and Eq. [A6] is linear and easily solved. APPENDIX
Ou \
The outward normal derivative on the boundary may be written
II
Terms herein retain their significance as previously defined. The first term in Eq. [5] and [A4] contains second-order differentials. These terms can be reduced to first-order differentials by the application of Green's lemma. If the domain {2 in the X - Y plane is bounded by a simple closed curve r that consists of a finite number of smooth arcs and if u(X, Y), v(X, Y), Ou/ OY, and Ov/OX are continuous at all points in f/and its boundary P then,
ff(
[B5]
Y (*re)
N{.ve)
dY
£
= fr (udX + vdY)dr.
[B1]
Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
FIG. B1. A section of the boundary r, showing the outward normal, N, and the angles 0x and Oy.
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION 0¢ _ 0¢ cos(0x) + 04~ cos(0r )
[B6]
ON OX -~ where cos Oxis the cosine of the angle between the positive X direction and the normal drawn outwards from F at the point where O¢/ON is computed; cos 0r has the corresponding meaning. Consider Fig. B 1 which shows the direction of the outward normal to the curve r at 0. d r is an infinitesimal length on r at 0, and is related to the infinitesimally small lengths dX and dY by the expression
d r 2 = dX 2 + dY 2.
dY dr
m
Equation [B11] is substituted into Eq. [5] which is written
(OPi Odp+ OPi O¢~axa Y Re
-ff(xVhChe-V°)dXdY ~e
/ m r \
+
F
, [B7]
When the derivatives dX/dF and dY/dP are computed with respect to the curve F they will have algebraic signs that depend upon the direction taken by F. For a line F surrounding ~ the convention is to take the forward F direction such that ~ is always on the lefthand side of F. Thus it follows that cos Ox
57
+ fr2Pi~NdF2=O.
[B12]
Recognizing that F is made up of r l and r 2 and that the integral along F1 is zero, then after multiplication by - 1 and simplification Eq. [B12] is written
[B81 9e
and cos 0y -
dX
dF"
[B9]
+ff,i~ ( ~ VhChe-Vh*)dXdY ~e
After multiplication by d r and substitution for cos Ox and cos Or Eq. [B6] is written as
o¢ d r = 0¢ d Y -
o-9
~
0¢
~ dX.
[B10]
Equation [B10] is substituted in Eq. [B5] which after some rearrangement may be written as
a2~ + ~-~) °2' ~!x,a"r f j "e['tox2 =
_ff t~_~_~_~+_~_~]lOPi 0¢ OPi O¢)dXd Y - fr Pi-d-~dr. 0¢
[Bll]
-fr2
e i ~[2 - - Nra ~ A i e\) d
F2 =
0.
[B13]
Computationally the evaluation of the line integration along F2 is straightforward, although the transformation between Cartesian and curvilinear coordinates is not readily apparent. In local curvilinear space a finite element will have nodes at each comer of a square with ~ and ~ coordinates of (-1, -1), (1, - 1), (1, l) and ( - 1, 1), respectively. For finite element adjacent to the boundary F2, the line integral in Eq. [B13] will need to be evaluated along one or more sides. For finite elements not adjacent to F2 the integral will be zero. The X and Y coordinates of a point on or within a finite element may be expressed in terms of the element shape funcJournal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
58
JAMES AND WILLIAMS
tions and nodal values element. That is
(Xi, Y;)
for that
X= ~P;Xi Y= ~ PiYi.
ing sides of the finite element are analogously defined as
[B14]
The chain rule of differentiation when applied to Eqs. [B14] yields m
m
~
Zi d~-~-
- ~ Zi an. [B15]
There will be two directions of integration in the transformed coordinates, one parallel to the ~ axis and the other parallel to the n axis. Thus it is possible to identify a particular line integral for each side of the square transformed finite element in ~-n space. First, there is the line integral along ~ = - 1 for = - 1 to ~ = + 1. In this situation the shape functions Pi(~, 7) are a function of ~ only, that is P,. = - 1 ) . This means that Eq. [BI5] is rewritten as
Pi(~,
dr=(XOP;Yg)d~
[B16]
+
along
n = +1
dr 2 =
along
~ -- - 1
A4d ~
[B19]
1. 2. 3. 4. 5. 6. 7. 8.
12.
[B17]
It is convenient to define this the boundary transformation along 71 = - 1 as
14.
[B18]
15. 16.
The transformations along the three remain-
17.
JournalofColloidandInterfaceScience,Vol. 107,No. 1, September1985
~ 1 Pi ~
1 Azd~ +
f+-iPi2--~AIe 1
A3d(
-~" f+-I1 ei ~ gO"
m4d~ [B20]
REFERENCES
13.
dP2 = Ald~.
d F 2 = A3d~
and remember that A~, A:, A3, or A4 will be nonzero only for sides of the finite element that are part of the boundary P2.
10. 11.
-~-X;
(i~OPi '271/2
2
9.
so that using Eq. [B7] we write dr2 =
~=-F1
where A2, A3, and A4 are determined in the same way as A1. We can now write the third term of Eq. [B13] as
dX=Q~X~)d~+Q~X~)d~I dY=
d F 2 = A2d~ along
Helmholtz, J. H., Wied. Ann. 7, 358 (1879). Stern, O., Z. Electrochem. 30, 508 (1924). Grahame, D. C., J. Chem. Phys. 18, 905 (1950). Gouy, G., J. Phys. 9, 457 (1910). Chapman, D. L., Philos. Mag. 25, 475 (1913). Nir, S., and Bentz, J., J. Colloid Interface Sci. 65, 399 (1978). Bentz, J., J. Colloid Interface Sci. 90, 127 (1982). Duchin, S. S., Derjaguin, B. V., and Semenikhin, N. M., Dokl. Akad. Nauk. SSSR 192, 357 (1970). Sigal, V. L., and Shama'nski, V. E., Dokl. Akad. Nauk. Ukr. SSR, B 4, 346 (1970). Stigter, D., J. Electroanal. Chem. 37, 61 (1972). Brenner, S. L., and Roberts, R. E., J. Phys. Chem. 77, 2367 (1973). Abraham-Schrauner, B., J. Colloid Interface Sci. 44, 79 (1973). White, L. R., J. Chem. Soc. Faraday II. 73, 577 (1977). MacGillivray, A. B., and Winkleman, J. J., J. Chem. Phys. 45, 2184 (1966). Mills, R. A., Biopolymers 9, 1511 (1970). Philip, J. R., and Wooding, R. A., J. Phys. Chem. 52, 953 (1970). Stigter, D., J. Colloid Interface Sci. 53, 296 (1973).
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION 18. van der Drift, W. P. J. T., de Kaiser, A., and Overbeek, J. Th. G., J. Colloid Interface Sci. 71, 67 (1979). 19. Torrie, G. N., and Valleau, J. P., Phys. Lett. 65, 343 (1979). 20. Debye, P., and Hiickel, E., Phys. Z. 24, 185 (1923). 21. Duke, G. P., Indian J. Phys. 17, 189 (1943). 22. Sparnaay, N. J., Rec. Tray. Chim. Pays-Bas 78, 680 (1959). 23. Hogg, R., Healy, T. W., and Fuerstenau, D. W., Trans. Faraday Soc. 62, 1638 (1966). 24. Muller, M., Phys. Z. 28, 324 (1927). 25. Guggenheim, E. A., Trans. Faraday Soc. 55, 1714 (1959). 26. Loeb, A. L., Overbeek, J. Th. G., and Wiersema, P. H., "The Electrical Double Layer Around a Spherical Colloidal Particle." MIT Press, Cambridge, Mass., 1961. 27. Hoskin, N. E., and Levine, S., Phil. Trans. Roy. Soc. A 248, 633 (1956). 28. McCartney, L. N., and Levine, S., J. Colloid Interface Sci. 30, 345 (1969). 29. Bresler, E., J. Colloid Interface Sci. 33, 278 (1970). 30. Desai, C. S., and Abel, J. F., "Introduction to the Finite Element Method." Reinhold, New York, 1972. 31. Galerkin, B. G., Vestn. Inzh. 19, 897 (1915).
59
32. Oden, J. T,, Adv. Astronautical Sci. 24, 17 (1967). 33. Oden, J. T,, Int. J. Num. Meth. Eng. 1, 247 (1969). 34. Szabo, B. A., and Lee, S. B., Int. J. Num. Meth. Eng. 1, 301 (1961). 35. Oden, J. T., "Finite Element Analysis of Nonlinear Continua." McGraw-Hill, New York, 1972. 36. Wylie, C. R., "Advanced Engineering Mathematics." McGraw-Hill, New York, 1960. 37. Zienkiewicz, O. C., "The Finite Element Method in Engineering Science." McGraw-Hill, London, 1971. 38. Irons, B. M., Int. J. Num. Meth. Eng. 2, 5 (1970). 39. Strang, G. S., and Fix, G. J., "An Analysis of the Finite Element Method." Prentice-Hall, Englewood Cliffs, N. J., 1973. 40. Brenner, S. L., and McQuarrie, D. A., J. Colloid Interface Sci. 44, 298 (1973). 41. James, A. E., and Williams, D. J. A., J. Colloid Interface Sci. 79, 33 (1981). 42. Flegmann, A. W., Goodwin, J. W., and Ottewill, R. H., Proc. Brit. Ceram. Soc. 13, 31 (1969). 43. Williams, D. J. A., and Williams, K. P., Trans. Z Brit. Cemm. Soc. 81, 98 (1982). 44. James, A. E., and Williams, D. J. A., Adv. Colloid Interface Sci. 17, 219 (1982). 45. Rall, L. B., "Computational solution of non-linear operator equations." Wiley, New York, 1962.
Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
Received March 8, 1984; acceptedJanuary 2, 1985 A Galerkin finite element scheme is used to provide flexible numerical solutions to the PoissonBoltzmann equation for electrical double layers in one and two dimensions. A Newton sequence is used for the solution of the set of nonlinear equations arising from the finite element discretization procedure. Convergence of the Newton sequence is rapid, generally occurring in less than four iterations. Isoparametricfinite elements are used for computational ease. Error analysis is applied to the solution of the linearized Poisson-Boltzmann equation (Debye-Hiickelapproximation) and also to the full iterative scheme. Numerical solutions are presented for (i) overlapping cylindrical and planar electrical double layers, and (ii) the overlapping electrical double layers of two interacting platelets. © 1985 Academic Press, Inc. INTRODUCTION
colloidal particles. Analytic solutions to the Poisson-Boltzmann equation are limited to special cases in one dimension. Complete solutions for the electrical potential decay from an infinite flat plate and have been given by Gouy (4) and Chapman (5) who consider the plate to be immersed in a binary symmetrical electrolyte. Nir and Bentz (6) and Bentz (7) have considered the effects of mixtures of electrolytes. The potential distribution for other geometries has been given for single spheres (Duchin et al. (8); Sigal and Shamanski (9); Stigter (10); Brenner and Roberts (11); Abraham-Schrauner (12) and White (13)) and single cylinders (MacGillivray and Winklemann (14); Duchin et al. (8); Sigal and Shamanski (9); Mills (15); Philip and Wooding (16); Stigter (17); van der Drift (18); Torrie and Valleau (19)). When the surface potential is 425 mV it is possible to linearize the exponential term in the Poisson-Boltzmann equation. This simplification is known as the Debye-Hiickel linearization and solutions to this equation are termed Debye-Hfickel solutions. The linearization removes the dependence of the solution on specific electrolyte types, so that Debye-Hiickel solutions are readily applica-
A knowledge of the distribution of electrical potential about colloidal particles dispersed in electrolyte is of great practical and theoretical importance. Typically, it is needed to calculate the free energy of interacting electrical double layers, which together with the attractive van der Waal's energy and energy attributable to hydrodynamic interactions, determine the stability of colloidal systems. The notion of the electrical double layer was proposed by Helmholtz (1), who likened it to a simple electrical capacitor. Electrical double layer theory has been modified by Stern (2) and Grahame (3). These authors split the electrical double layer into two regions: the inner or fixed region and the outer or diffuse region. The former is determined by the closest approach of counterions and the variation of potential in the region is likened to that in an electrical capacitor, In the diffuse layer, the distribution of potential is described by the Poisson-Boltzmann equation. The Poisson-Boltzmann equation is the fundamental differential equation describing the distribution of electrical potential about 44 0021-9797/85 $3.00 Copyright © 1985 by Academic Press, Inc. All rights of reproduction in any form reserved.
Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION
ble to electrolyte mixtures. Debye and Hiickel (20) give a solution for a single spherical electrical double layer, while Dube (21) and Sparnaay (22) give results for a single cylindrical double layer and two overlapping cylindrical double layers, respectively. The case of two parallel plates has been analysed by Hogg, Healy, and Fuerstenau (23). Numerical solutions to the Poisson-Boltzmann equation have been given for a single sphere (Muller (24); Guggenheim (25), and Loeb, Overbeek and Wiersema (26)), two identical spheres in a binary symmetrical electrolyte (Hoskin (27); McCartney and Levine (28)), a single cylinder (Stigter (17)), and the potential distribution between two parallel plates in mixed ionic systems (Bresler (29)). Analytic solutions to the Poisson-Boltzmann equation are limited to regularly shaped individual colloidal particles. Similar restrictions also apply to existing numerical solutions. Numerical solutions of the PoissonBoltzmann equation have used the finite difference method combined with first-order iterative methods and have not been flexible enough to allow for changes in size, shape, and geometry of the system under consideration. Currently most analytic and numerical solutions to the Poisson-Boltzmann equation have been one-dimensional simplifications which have necessarily limited application in real situations. There is need for a general numerical solution to the Poisson-Boltzmann equation for arbitrary geometry. To this end, a finite element scheme is presented using a second-order iterative scheme to speed the convergence of the nonlinear set of equations obtained from the discretization process. 2. THEORY
The Poisson-Boltzmann equation written for a two-dimensional electrical double layer is 020
020
~ 2 + 0 Y~ = -0(0) OX
[ 1]
where X and Y are nondimensional distances,
45
0 is the nondimensional electrical potential, and p(O) is the nondimensional space charge density. That is X = rx, Y = ry, and ¢ = op/ kT, where x and y are Cartesian coordinate distances, K = (2NAIe2/~kT) 1/2, e = electronic charge, ~b = electrical double layer potential, k = the Boltzmann constant, T = absolute temperature, ~ = permittivity of the electrolyte, NA = Avogadro's number, I = (ZV2hCh)/2 = the ionic strength of the electrolyte, with vh and ch being the valency and concentration of the hth ionic species. The summation is carried out over all ionic species in the electrolyte and vh contains the sign of the charge of the hth ionic species. Consistent S.I. units have been used. For example the concentrations ch and I are in kmole m -3 and K the DebyeHfickel parameter has units of m -1. The nondimensional space charge density is written
p(O) = (r--g---ekT) ~ VheNAche-"h* 1
- - - Z VhChe-vh*. 21
[2]
The summation again is taken over all ionic species present in the electrolyte. In order to solve Eq. [ 1] it is necessary to define the potential, 0, and the potential gradient 00/0N at the boundaries, 1", of the chosen solution domain, ft. Here N (=Kn) is the outward normal direction from the boundary 1". The boundary condition is written as ¢ = 0b on r I [3] and 00
--KO"
ON - 2NAIe on 1"2
[4]
where F1 includes all parts of F where the nondimensional potential is fixed, 1,2 includes all parts of F where the nondimensional potential gradient is defined, and a is the surface charge density (C m -E) on 1"2, (ga/ 2NAIe) is the nondimensional surface charge density. In order to produce a numerical solution to the Poisson-Boltzmann equation using Journal of Colloid and Interface Science,
Vol. 107,No. 1, September1985
46
JAMES A N D W I L L I A M S
the finite element technique it is necessary to divide the solution domain into a grid (as in the finite difference procedure). This grid need not be regular and the finite elements thus produced are usually triangular or rectangular. In the present work the latter type were used. Each finite element has a fixed number of nodal points along its boundary and the numerical solution is generated at these points. The variation of a parameter (e.g., 4, X, Y) within a finite element is described in terms of the nodal values. Thus for a given area, the accuracy of the approximation increases as the number of finite elements increases and their size decreases. The finite element method is an extension of the Rayleigh-Ritz-Galerkin technique, whereby the functional form of the original differential equation is discretized by assuming the solution to be given by sets of overlapping polynomials whose nodal values have to be determined. The discretization process is applied to each finite element in the solution domain so that the variation of the unknown variable in each finite element is described by a set of polynomials which are known as shape functions. The order of the polynomial used characterizes the type of finite element by specifying the number of nodal points on the side of a finite element. The discretized set of equations describing the variation of the unknown parameters in individual finite elements is assembled using the direct stiffness method (Desai and Abel (30)) into an overall or global set, which is subsequently solved for the nodal values of the unknown variables. In the present work the Galerkin (31) method is used to generate the equivalent functional form of the Poisson-Boltzmann equation, and this is subsequently discretized. The usefulness of this method was recognized by Oden (32, 33) and Szabo and Lee (34); further, Oden (35) suggests that it is the most powerful technique for generating finite element representations of nonlinear differential equations. The Galerkin method does not operate directly with the relevant differential Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
equation but with the equivalent functional form of the continuous problem, of which the solution is assumed to be a combination of trial functions, say Pi- Thus at any node i within a finite element having a total of m nodes it follows that the equivalent functional form of Eqs. [1], [3], and [4] is y,/d2q5
1
02~b
VhChe_Vh+)
f~e
× d~2e + f~ Pi(X, Y)(c~ -
4)b)dF~
dr
I Ka
04,]dF~
wherein Pi(X, Y) is the weighting function with respect to the ith node and the superscript e refers to a single finite element. Green's lemma (see Wylie (36)) is used to convert the second-order partial derivatives in Eq. [5] into products of first-order terms (see Appendix II for details). The transformed equation is written
~e
,ax o x
& 2I (~ VhChe-v*~)dXdg Kff
e
+ ~r Pi(2-~Aie)dr2=O.
[6]
It is worthwhile noting that the second term in Eq. [5] vanishes as (4) - 4)b = 0) by definition on r~, while the third term in Eq. [6] is the line integral of the nondimensional surface charge along F~. The first two terms are evaluated over the area of the finite element and implicitly contain boundary conditions of the form 4~ = q~b since this condition is specified at any node on the boundary. The final step in the finite element formulation is the discretization of the equivalent functional form of the Poisson-Boltzmann equation given by Eq. [6]. Within any
47
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION
particular finite element the variables X, Y, and ¢ can be written in terms of their nodal values as x e : E pe(~, n)X e = [e]e{x}e
ye = X pe(~, y)y~
=
[p]e{y}e
q~e = E pe(~, n)¢e = [p]e{¢}e.
[7]
charge, % along each side of a rectangular finite element in turn. This is consistent with the convention that the solution domain is always to the left of the line along which the integral is performed. The values A, (i = 1, 2, 3, 4), are given by Eqs. [B17]-[B19] while the Jacobian determinant IJI --- O(S, Y)/O(~, n). These parameters are required to change from the cartesian (X, Y) coordinate system to the curvilinear (~, n) coordinate system in which the integration of the equivalent functional form of the Poisson-Boltzmann equation is performed. Equation [8] can conveniently be incorporated into a matrix equation for the m nodes in a finite element [F]e{dp}e+ {C(~b)}e-[ - { S } e = { 0 } e [9]
The summations are made over the m nodes in the finite elements; [p]e is the (1 × m) line vector of shape functions; {X}, {Y}, and {¢} are the (m × 1) column vectors of the nodal values of X, ]1, and ¢, respectively. P~((, n) is a polynomial of chosen degree (usually less than 3) and ( and n are curvilinear coordinates introduced to enable the integration to be performed over a uniform region. If the same polynomials, Pg, are used where [F]e is an (m × m) symmetric positive for both the weighting function in Eq. [5] definite matrix of which the general member and also for the description of a variable f0 is written within a finite element then the finite elements are isoparametric (Zienkiewicz (37)). f/j = :__'-11__'-I 1 ~ 1 Equation [7] is substituted into Eq. [6] to give (OPiOPj OPiOPJ\jd d
×\OXOX+-~-~)l
f~'-I f~"l (Opi O[p]e Opio[p]e ) 1 l ~OX OX {l#}e"~ OY OY {~)}e
[101
The (m × 1) column vector {C(¢)} contains the nonlinear terms related to the space charge density and has the form at the ith row of
f~'-lf-lpi 1 I 2I
× lJld~dn -
[ ~ 7/.
× (X VhChexp(-- Vh[Ple{¢}e))lJldtidrl cR¢)
= 1
-[-
Pi E (eql~ ql)mld~ 1 ql
+
£1Pi(~ Pq2~q2)A2dn
+
j
1
1
× ( Z VhChexp(--Vh[p]e{¢}e))lJld~dvl [1 1] and the surface charge boundary condition vector, (m × 1) has the general term
q2
~- I +1
L
Se=
Pi(Z Pq3~q3)A3d~
Pi(~_,Pql~ql)Ald~
I
q3
ql
j
"+ 1
+ +
Pi(E Pe---e)A4dn +1 q4
= 0
[81
where ~q = (K/2NAIe)~q, and ---q = 0 if q is not a boundary condition node. The summations ~q~, ~]q2, ~q3, and Zq4 refer to the evaluation of the nondimensional surface
-1
+ J+l
Pi( ~ Pq2==-q2)~2d~ q2
Pi(~ eq3~q3)m3d~ q3
+
f? P,(Z Pe%4),~dn.
[12]
1 q4 Journalof ColloidandInterfaceScience,Vol. 107, No. 1, September1985
48
JAMES AND WILLIAMS
This last term assumes nonzero values only at nodes where surface charge is specified as a boundary condition. Because of the nonlinearity of terms in the space charge vector Eq. [9] cannot be inverted immediately. In the present work, a Newton sequence (see Rail (45)) is used to solve the nonlinear set of equations described by Eq. [9]. Let a nonlinear operator A~(~) be defined as
A~(4~) = [Fle{~b} + {C(q~)} eq-
{S} e
[131
It is apparent from Eq. [18] that [ J A G ] e is symmetrical, which is useful in the numerical computations. The actual Newton sequence required is now written as
[SAC]e{ Adp}e*+l = - [ F y { q ~ } e . - { C ( 4 , ) } e . - {S}~n. {4~}~n.+1 = {~}e. + {dX4~}e.+l
[19]
3. COMPUTATIONAL ASPECTS OF THE SOLUTION
SO that the Jacobian matrix A~(q~) is written The Poisson-Boltzmann equation has been discretized and a suitable Newton sequence OA~ck Ael(dP) - - - [ J A C ] e. [141 has been presented for the solution of the resultant set of nonlinear equations describing the potential distribution within a single finite The Newton sequence required to solve element. The solution for the potential A~(qS) = 0 is defined by throughout an entire region which has been A~(~b){A~b}e*+, = -A~(4~),* [15] divided into finite elements is obtained by constructing global matrices corresponding {q~}e.+l = {q~}ne. q- {Aq~e}n,+l to [JAC] e and [F] e and the global column where the subscript n* refers to values of ~b vectors {C(q~)}e and { S } e. This global assemand 2x~ at the n*th iteration. An initial set bly is achieved using the direct stiffness of values to start the sequence is taken from method (e.g., see Oden (35) or Desai and the solution to the Debye-Hiickel problem Abel (30)). A frontal solution scheme (see (see Appendix I). The convergence of Newton Irons (38)) is used to invert the global set of sequences is discussed in detail by Rail (38). equations. The Newton sequence does not Practically, A~(q~) is an m x m matrix, the require the recalculation of [F] e on each iteration. The coefficients in every other mageneral member being given by trix and vector have to be evaluated at each O(a~(ck)) jac~ - - [ 16] step. Provided a good initial guess is made 0~j using the Debye-Hfickel solution (see ApHere a~i(ck) is the i th row of A~(q~) such that, pendix I) the Newton sequence usually converges within three iterations. Programs imusing Eq. [ 11 ], plementing this numerical scheme have been m written in F O R T R A N for an ICL 1904S aoi = ~ f,~q~j + ce(qs) + s e. [17] computer. as
j=l
Thus, using Eqs. [10], [11], [16], and [t7] m and noting that q5 = ~i=l PigJi, 1 f+l
jac~=f°
e__
~d-i
i'+l
,-i
pipj
p*
x ( Z VhChexp(-vhck))lJId~dn.
[18]
h=l Journal of Colloid and Interface Science, Vol. 107,No. 1, September1985
4. ERROR ANALYSIS
There are no analytic solutions to the Poisson-Boltzmann equation in two dimensions even for such a simple geometry as a single isolated charged cylinder in a binary symmetrical electrolyte. Although this is a one-dimensional problem in cylindrical coordinates, it would have provided a suitable
NUMERICAL
SOLUTION
OF THE POISSON-BOLTZMANN
test in rectangular Cartesian coordinates. There are a number of approximate analytic solutions but these are of little use in ascertaining the accuracy of the finite element solution. In one dimension there is an analytic solution to the Poisson-Boltzmann equation (Gouy (4) and Chapman (5)) for binary symmetrical electrolytes. The equation describing the decay of potential from an infinite fiat plate into an infinite fiat double layer is written as e~°/2 + 1 + (e ~°/2- l)e - x 4 ' = l n e ~ o / 2 + 1 - ( e *°/2- l)e - x "
[20]
Here 4)0 is the reduced surface potential. For 4'0 < 1 the corresponding Debye-Hfickel solution is [21 ]
4) = 4'oe-x.
Equations [20] and [21] were derived using the boundary conditions X=0,
4'=
4,0
X= ~,
~-~=0.
a4'
[22]
The infinite boundary condition is problematical in some numerical solutions. In the present work, it was sufficient to choose a value of X large enough that d4'/dX _~ O. In practice, this occurs at around 30 units from any surface. The boundary conditions used in the flat plate solutions are X=O
4'=4'0
X = 30,
d4'/dX = 0.
49
n is the order of the functional that has been constructed from the original differential equation, and B is a constant. In the present work n = 1, and k = 2, 3, 4, corresponding to linear, quadratic, and cubic finite elements. Equation [24] suggests that the error in each of these finite elements is proportional to h 2, h 4, and h 6, respectively. The errors between the finite element solution to the distribution of electrical potential in a plane double layer and Eqs. [20] and [21] have been computed and are shown in Figs. 1 and 2. The values of B, obtained by regression of ]error[ o n h (2K-n) are shown in Table I. Although the finite element solution to the Poisson-Boltzmann equation uses an iterative scheme the values of B are not significantly different from those seen for the linear Debye-Hfickel solution. It is notable that the Newton sequence takes no more than three iterations to converge while other iterative procedures used to solve the PoissonBoltzmann equation have used over 20 iterations (see Hoskin and Levine (27); McCartney and Levine (28)). 5. E R R O R S
IN TWO-DIMENSIONAL SOLUTIONS
There are no suitable analytic solutions to test the convergence of the numerical solution 0-
-2-
[231
The forcing of the potential gradient at large X does not introduce any significant errors into the computer solution especially as an iterative technique is used. Strang and Fix (39) have shown that the simplest and most fundamental error estimate to linear finite element problems is given by lerrorl --< B2h 2(k-n)
EQUATION
[241
where h is the length of the element side, k is the order of the shape functions (P;((, n)),
L
o
_/,-
-6
-8
-10 -2 log
(h)
FIG. 1. The relationshipbetween [error[and h, length of a finiteelement,when linear(LFES),quadratic (QFES), and cubic(CFES)finiteelementsare used in the numerical solution of the Poisson-BoltzmannEquation. Journal of Colloid and Inte(face Science, Vol. 107, No. 1, September 1985
50
JAMES AND WILLIAMS 0
Figure (3) shows the arrangement of finite elements and relevant boundary conditions -2 for the numerical solution to the PoissonBoltzmann equation about a single cylinder. Typical results for the Debye-HiJckel solution (Eq. [24]) and the numerical solution using -6 quadratic finite elements, are shown in Table o= II. It can be seen in Fig. 3 that for the aspect -8 ratios of the elements to be kept in the region of unity, it is necessary to increase the density of elements (i.e., reduce their size) as the -10 -2 cylinder wall is approached. In practice, it log (h) was found that it was not necessary for the FIG. 2. The relationship between [errort and h, the aspect ratio of elements close to the cylinder length of a finite element, when linear (LFES), quadratic to be close to unity for good numerical (QFES), and cubic (CFES) finite elements are used in the numerical solution of the Poisson-Boltzmann Equa- solutions to be obtained. The main criteria for the simple single cylinder geometry seems tion. to be the radial length of the element, which once small enough, ensures good agreement to the Poisson-Boltzmann equation with re- between analytic (Eq. [25]) and numerical spect to element size for a two-dimensional schemes. This result is possibly linked to the problem. It would seem reasonable that the constraint that for a single cylinder &b/ON relationship between the length of an element = 0, where N is the direction normal to the side and the absolute error will hold in both radial direction. When the geometry becomes Cartesian (X, Y) directions if finite elements more complex, e.g., two or more cylinders, with an aspect ratio of unity are used (ratio Odp/ON d= O. of m a x i m u m length to m i n i m u m length, see Desai and Abel (30)). A major attraction of 6. TWO-DIMENSIONAL NUMERICAL the finite element method is its flexibility in SOLUTIONS allowing elements of different size and shape One-dimensional solutions to the Poissonto be used in the same mesh; clearly, this attraction is lost if all the elements are equi- Boltzmann equation have limited applications sized. Further, in m a n y two-dimensional in colloid chemistry. In m a n y instances it is problems, it is not possible to use rectangular necessary to know the electrical potential elements with an aspect ratio of unity. It can be shown (Sparnaay (22)) that under the restrictive conditions of the D e b y e TABLE I H/ickel approximation, an analytic solution Values of B, Determined by Linear Regression of the can be obtained. That is Error Relationship (Eq. [25] and Fig. 1) for Linear, Qua4~ = 4~oKo(R)/Ko(Ro) [25] dratic, and Cubic Finite Elements
j
where R is the nondimensional radial distance from the axis of the cylinder and R0 is the cylinder radius. ~bo is the reduced surface potential of the cylinder. K0( ) is a modified Bessel function of the second kind and zero order. Journal of Colloid and Interface Science, Vol. 107, No. I, September 1985
Linear element Quadratic element Cubic elements
Debye-Hfickel solution
PoissonBoltzmann solution
5.46 × 10-2 7.03 × 10-3 10-3
7.16 X 10 -2 7.41 X 10-3 3.16 × 10-4
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION
element
51
solution
FIG. 3. Arrangement of finite elements and relevant boundary conditions for a single cylinder.
distribution in the double layers about two or m o r e particles immersed in electrolyte. Some interesting and relevant geometries are those involving cylindrical or plate-like particles. A two-dimensional solution to the P o i s s o n - B o l t z m a n n equation can be utilized in such cases and others where s y m m e t r y can be used to reduce the dimensionality o f the problem. The problem o f interacting particles having different geometries m a y either be a real case where the geometry is known, e.g., between pairs o f tobacco mosaic virus (cylinder-cylinder) (see Brenner and McQuarrie (40); James and Williams (41)), between protein molecules and a cell wall (cylinder-plate), or in a less well defined system, where some geometric shape m u s t be assumed to evaluate interparticle electrical double-layer potentials,
e.g., kaolinite suspensions (Flegmann, G o o d win, and Ottewill (42); Williams and Williams (43); James and Williams (44)). Firstly, the problem o f overlapping cylindrical and planar electrical double layers typifying the protein cell wall interaction is examined. Figure 4 shows the b o u n d a r y condition o f the cylinder plate geometry. The general layout o f the finite element mesh is shown in Fig. 5. A n u m b e r o f mesh refinements were required to obtain a convergent solution. The areas o f modification are shown in Fig. 5. Figures 6 and 7 show results for the electrical double-layer interaction between a cylinder and a plate at different separations. The increase in potential gradient reflecting an increase in the attractive free energy is apparent at the smaller separation. The distortion o f the potential distribution about the
TABLE II Comparison between Debye-Hiickel Solution for a Cylinder and the Analytic Debye-Hiickel Solution (Eq. [25]) for the Potential about a Cylinder, Using Quadratic Finite Elements~ Number of elements R
26
39
49
56
66
Anflytic equation [25]
0.2 0.7 1.2 2.2 4.2 6.2
2.0 0.806 0.380 0.107 1.87 X 10 -2 2.07 X 10 -3
2.0 0.759 0.367 0.103 1.02 X 10 -2 1.15 X 10-3
2.0 0.755 0.364 0.102 1.02 X 10-2 1.15 X 10-3
2.0 0.7543 0.3637 0.1019 1.019 X 10 -2 1.145 X 10-3
2.0 0.7538 0.3634 0.1019 1.019 X 10 -2 1.141X 10-3
2.0 0.7538 0.3634 0.1019 1.019 X 10 2 1.141X 10 -3
a
Cylinder radius = R0 = 0.2. Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
52
JAMES
AND
9¢ - - =0 9N
to describe the interaction and geometry of two interacting platelets, which could be clay particles. Typical solutions are illustrated in Figs. 10 and 11. The potential gradient is seen to increase at the closer separation (Fig. 11). It is also clear that the electrical field associated with the edge becomes more localized and distorted at close separation (cf. Fig. 9). Both of these observations are consistent with an increase in the attractive free energy of the system.
y=÷-~
¢
/
9~ =0 9N X= - e e
:
WILLIAMS
~02
¢=~0
7. S U M M A R Y
9 ~ = 0 along AA' 9N /
S /
9~ =0 9N
--
V = -~ FIG. 4. B o u n d a r y
conditions for the cylinder-plate
geometry.
cylinder shows the effects of an increasingly strong interaction at close separation. Figures 8 and 9 show the boundary conditions and basic finite element mesh used
;
This work has shown that in one or two dimensions, the finite element method is ideally suited to solutions of the nonlinear Poisson-Boltzmann equation. The variable geometry of interacting charged surfaces is more readily accommodated by a flexible, nonuniform finite element mesh than by a finite difference scheme. The finite element procedure produces a set of nonlinear equations which are solved in the present work with the aid of a Newton sequence; this sequence converges rapidly when a trial so-
~N =0
3
3
3
3
3
3
o2
1
First
modification 2
Second modification
3
Third modificotion
+o
9N =0
Ol
--=0 9N
FIG. 5. Finite e l e m e n t m e s h for cylinder-plate geometry.
Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
53
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION
-5.0 -L.O
-30
-2.0
-1.0
0.0
1.0
2.0
3.0
4.0
5,0
6.0
7.0
8.0
9.0
X
FIG. 6. Distribution of potential ~b, for the interaction of a cylinder (Ro = 1.0) with a plate at a separation X = 8.0.
lution obtained from the linearized DebyeHiickel equation is used. Strang and Fix (39) have given a criterion for the convergence of finite elements in linear problems (e.g., the Debye-Hiickel equation). Interestingly the present work has shown that the error relationship is applicable to the nonlinear Poisson-Boltzmann equation, when either linear, quadratic, or cubic finite elements are used. The finite element method is a very convenient and powerful method of studying the various modes of potential interaction that can occur between particles of irregular geometry. In contrast to finite difference schemes the actual geometrical surface of the particles can be easily incorporated in the numerical scheme. The technique is readily extended to three-dimensional problems such
as the potential interactions between pairs of spheres and to more complicated systems, such as arrays of particles. APPENDIX I: FINITE ELEMENT FORMULATION FOR THE DEBYE-HOCKEL APPROXIMATION
The Debye-Hiickel approximation to the two-dimensional Poisson-Boltzmann equation is written 0%
0% + ~ = 4,
[A1]
with boundary conditions ~b = 4~b on P L
Odp ON
[A21
K~
- -
2NAIe
o n 172.
[A3]
Journal of Colloid and Interface Science, Vol. 107, No. 1, Septembert985
54
JAMES A N D WILLIAMS
7.0 6.0 5.0 &0
Note that this equation differs from Eq. [6] in that the second term contains no nonlinear function of 4. Substituting for 4~ and using isoparametric finite elements leads to the matrix equation for a single finite element
3.0
[Dle{dp}e + {S} e =
0,
[A6]
20 1.0 0.0 -19 -2.0
where
+l (Op~opj OPiOPj ) + -if-f-if-f+ PiP/ 4J= f-l f-, tOX OX x IJld~n [A7] +1
-3.0 -/,.0 -5.0
and
Si =
£.-1Pi(Z Pql-~ql)Ald~ +F' 1
-6,0 -70
ql
1 -5.0
-40
-3.0 -2.0
-1.0
00
1,0
20
30
4.0
q2
5.0
X
Pi(
FIG. 7. Distribution of potential 4, for the interaction of a cylinder (R0 = 1.0) with a plate at a separation X = 4.0.
1
Z
Pq3m-q3)A3d~
q3
j
,t_ 1
+
Pi(~ Pq4=--q4)A4dT1• [A8] +1
q4
Galerkin's method is applied to Eq. [A1][A3], so that at any node i in a finite element, it follows that / / / / /
Y =~
f f Pi,ox2 0y2-
9N
¢=~
X=0
+ fr~ P~(4a- 4,b)dr[
~N
ra're 2 = 0. [A4]
+fr~Pi
[oei 04
% 0/
¢~-,~
A'
/
+ ~ ) d
After applying Green's lemma and some rearrangement (see Appendix II), Eq. [A4] is written
ff
A
/ / / /
ori~y)dfU+ ff e''d"e
£e
£e
+
dr2 = 0.
[AS]
Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
/
¢ / / / / /
X = 0 9~
:0
~N
¢ 9~=0 ~N
FIG. 8. Boundary conditions for the interaction of two plates in edge-face mode.
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION
9~
9N
55
-0 02
02
01
~ "~N
0
FIG. 9. Finite element mesh for two interacting plates in edge-face mode.
4.0
5.0
6.0
7.0
8.0
9.0 10.0 11.0 12~ 13.0 l&0 15.0 16.0 1Z0
X FIG. 10. D i s t r i b u t i o n o f potential th, f o r the edge-face interaction o f t w o plates at a separation X = 8.0.
Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
56
JAMES
0.02 0.05 0:1 02
AND
WILLIAMS
Proof of Eq. [B1] may be found in most mathematical texts concerning line, surface, and volume integrals. (For example, see Wylie (36) p. 568 et seq.) On substitution of
05 1C ZO 50 5.0 4.0 3.0
v =e/0-0q~
20 1.0
T
0.0 -1.0
Y
[B21
and
04
[B3]
u = - P i O---Y Eq. [B 1] is written
-20 -3.0 -4.0
£
-5.0
o
o2
r
= fr ( - e i - ~ dX + Pz ox dY) d ,
-6.0
[B41
-ZO
so that 3.0
&0
50
6.0 X
7.0
&0
02, oei0,+
g.0 10.0 11.0 12.0 13.0
=,
F I G . 11. Distribution of potential ~, for the edge-face interaction of two plates at a separation X = 4.0.
£
04
04
× dXdY = fr Pi(ff-x dY - -~ dX ) • The matrix [D] e is symmetrical and Eq. [A6] is linear and easily solved. APPENDIX
Ou \
The outward normal derivative on the boundary may be written
II
Terms herein retain their significance as previously defined. The first term in Eq. [5] and [A4] contains second-order differentials. These terms can be reduced to first-order differentials by the application of Green's lemma. If the domain {2 in the X - Y plane is bounded by a simple closed curve r that consists of a finite number of smooth arcs and if u(X, Y), v(X, Y), Ou/ OY, and Ov/OX are continuous at all points in f/and its boundary P then,
ff(
[B5]
Y (*re)
N{.ve)
dY
£
= fr (udX + vdY)dr.
[B1]
Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
FIG. B1. A section of the boundary r, showing the outward normal, N, and the angles 0x and Oy.
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION 0¢ _ 0¢ cos(0x) + 04~ cos(0r )
[B6]
ON OX -~ where cos Oxis the cosine of the angle between the positive X direction and the normal drawn outwards from F at the point where O¢/ON is computed; cos 0r has the corresponding meaning. Consider Fig. B 1 which shows the direction of the outward normal to the curve r at 0. d r is an infinitesimal length on r at 0, and is related to the infinitesimally small lengths dX and dY by the expression
d r 2 = dX 2 + dY 2.
dY dr
m
Equation [B11] is substituted into Eq. [5] which is written
(OPi Odp+ OPi O¢~axa Y Re
-ff(xVhChe-V°)dXdY ~e
/ m r \
+
F
, [B7]
When the derivatives dX/dF and dY/dP are computed with respect to the curve F they will have algebraic signs that depend upon the direction taken by F. For a line F surrounding ~ the convention is to take the forward F direction such that ~ is always on the lefthand side of F. Thus it follows that cos Ox
57
+ fr2Pi~NdF2=O.
[B12]
Recognizing that F is made up of r l and r 2 and that the integral along F1 is zero, then after multiplication by - 1 and simplification Eq. [B12] is written
[B81 9e
and cos 0y -
dX
dF"
[B9]
+ff,i~ ( ~ VhChe-Vh*)dXdY ~e
After multiplication by d r and substitution for cos Ox and cos Or Eq. [B6] is written as
o¢ d r = 0¢ d Y -
o-9
~
0¢
~ dX.
[B10]
Equation [B10] is substituted in Eq. [B5] which after some rearrangement may be written as
a2~ + ~-~) °2' ~!x,a"r f j "e['tox2 =
_ff t~_~_~_~+_~_~]lOPi 0¢ OPi O¢)dXd Y - fr Pi-d-~dr. 0¢
[Bll]
-fr2
e i ~[2 - - Nra ~ A i e\) d
F2 =
0.
[B13]
Computationally the evaluation of the line integration along F2 is straightforward, although the transformation between Cartesian and curvilinear coordinates is not readily apparent. In local curvilinear space a finite element will have nodes at each comer of a square with ~ and ~ coordinates of (-1, -1), (1, - 1), (1, l) and ( - 1, 1), respectively. For finite element adjacent to the boundary F2, the line integral in Eq. [B13] will need to be evaluated along one or more sides. For finite elements not adjacent to F2 the integral will be zero. The X and Y coordinates of a point on or within a finite element may be expressed in terms of the element shape funcJournal of Colloid and Interface Science, Vol. 107, No. 1, September 1985
58
JAMES AND WILLIAMS
tions and nodal values element. That is
(Xi, Y;)
for that
X= ~P;Xi Y= ~ PiYi.
ing sides of the finite element are analogously defined as
[B14]
The chain rule of differentiation when applied to Eqs. [B14] yields m
m
~
Zi d~-~-
- ~ Zi an. [B15]
There will be two directions of integration in the transformed coordinates, one parallel to the ~ axis and the other parallel to the n axis. Thus it is possible to identify a particular line integral for each side of the square transformed finite element in ~-n space. First, there is the line integral along ~ = - 1 for = - 1 to ~ = + 1. In this situation the shape functions Pi(~, 7) are a function of ~ only, that is P,. = - 1 ) . This means that Eq. [BI5] is rewritten as
Pi(~,
dr=(XOP;Yg)d~
[B16]
+
along
n = +1
dr 2 =
along
~ -- - 1
A4d ~
[B19]
1. 2. 3. 4. 5. 6. 7. 8.
12.
[B17]
It is convenient to define this the boundary transformation along 71 = - 1 as
14.
[B18]
15. 16.
The transformations along the three remain-
17.
JournalofColloidandInterfaceScience,Vol. 107,No. 1, September1985
~ 1 Pi ~
1 Azd~ +
f+-iPi2--~AIe 1
A3d(
-~" f+-I1 ei ~ gO"
m4d~ [B20]
REFERENCES
13.
dP2 = Ald~.
d F 2 = A3d~
and remember that A~, A:, A3, or A4 will be nonzero only for sides of the finite element that are part of the boundary P2.
10. 11.
-~-X;
(i~OPi '271/2
2
9.
so that using Eq. [B7] we write dr2 =
~=-F1
where A2, A3, and A4 are determined in the same way as A1. We can now write the third term of Eq. [B13] as
dX=Q~X~)d~+Q~X~)d~I dY=
d F 2 = A2d~ along
Helmholtz, J. H., Wied. Ann. 7, 358 (1879). Stern, O., Z. Electrochem. 30, 508 (1924). Grahame, D. C., J. Chem. Phys. 18, 905 (1950). Gouy, G., J. Phys. 9, 457 (1910). Chapman, D. L., Philos. Mag. 25, 475 (1913). Nir, S., and Bentz, J., J. Colloid Interface Sci. 65, 399 (1978). Bentz, J., J. Colloid Interface Sci. 90, 127 (1982). Duchin, S. S., Derjaguin, B. V., and Semenikhin, N. M., Dokl. Akad. Nauk. SSSR 192, 357 (1970). Sigal, V. L., and Shama'nski, V. E., Dokl. Akad. Nauk. Ukr. SSR, B 4, 346 (1970). Stigter, D., J. Electroanal. Chem. 37, 61 (1972). Brenner, S. L., and Roberts, R. E., J. Phys. Chem. 77, 2367 (1973). Abraham-Schrauner, B., J. Colloid Interface Sci. 44, 79 (1973). White, L. R., J. Chem. Soc. Faraday II. 73, 577 (1977). MacGillivray, A. B., and Winkleman, J. J., J. Chem. Phys. 45, 2184 (1966). Mills, R. A., Biopolymers 9, 1511 (1970). Philip, J. R., and Wooding, R. A., J. Phys. Chem. 52, 953 (1970). Stigter, D., J. Colloid Interface Sci. 53, 296 (1973).
NUMERICAL SOLUTION OF THE POISSON-BOLTZMANN EQUATION 18. van der Drift, W. P. J. T., de Kaiser, A., and Overbeek, J. Th. G., J. Colloid Interface Sci. 71, 67 (1979). 19. Torrie, G. N., and Valleau, J. P., Phys. Lett. 65, 343 (1979). 20. Debye, P., and Hiickel, E., Phys. Z. 24, 185 (1923). 21. Duke, G. P., Indian J. Phys. 17, 189 (1943). 22. Sparnaay, N. J., Rec. Tray. Chim. Pays-Bas 78, 680 (1959). 23. Hogg, R., Healy, T. W., and Fuerstenau, D. W., Trans. Faraday Soc. 62, 1638 (1966). 24. Muller, M., Phys. Z. 28, 324 (1927). 25. Guggenheim, E. A., Trans. Faraday Soc. 55, 1714 (1959). 26. Loeb, A. L., Overbeek, J. Th. G., and Wiersema, P. H., "The Electrical Double Layer Around a Spherical Colloidal Particle." MIT Press, Cambridge, Mass., 1961. 27. Hoskin, N. E., and Levine, S., Phil. Trans. Roy. Soc. A 248, 633 (1956). 28. McCartney, L. N., and Levine, S., J. Colloid Interface Sci. 30, 345 (1969). 29. Bresler, E., J. Colloid Interface Sci. 33, 278 (1970). 30. Desai, C. S., and Abel, J. F., "Introduction to the Finite Element Method." Reinhold, New York, 1972. 31. Galerkin, B. G., Vestn. Inzh. 19, 897 (1915).
59
32. Oden, J. T,, Adv. Astronautical Sci. 24, 17 (1967). 33. Oden, J. T,, Int. J. Num. Meth. Eng. 1, 247 (1969). 34. Szabo, B. A., and Lee, S. B., Int. J. Num. Meth. Eng. 1, 301 (1961). 35. Oden, J. T., "Finite Element Analysis of Nonlinear Continua." McGraw-Hill, New York, 1972. 36. Wylie, C. R., "Advanced Engineering Mathematics." McGraw-Hill, New York, 1960. 37. Zienkiewicz, O. C., "The Finite Element Method in Engineering Science." McGraw-Hill, London, 1971. 38. Irons, B. M., Int. J. Num. Meth. Eng. 2, 5 (1970). 39. Strang, G. S., and Fix, G. J., "An Analysis of the Finite Element Method." Prentice-Hall, Englewood Cliffs, N. J., 1973. 40. Brenner, S. L., and McQuarrie, D. A., J. Colloid Interface Sci. 44, 298 (1973). 41. James, A. E., and Williams, D. J. A., J. Colloid Interface Sci. 79, 33 (1981). 42. Flegmann, A. W., Goodwin, J. W., and Ottewill, R. H., Proc. Brit. Ceram. Soc. 13, 31 (1969). 43. Williams, D. J. A., and Williams, K. P., Trans. Z Brit. Cemm. Soc. 81, 98 (1982). 44. James, A. E., and Williams, D. J. A., Adv. Colloid Interface Sci. 17, 219 (1982). 45. Rall, L. B., "Computational solution of non-linear operator equations." Wiley, New York, 1962.
Journal of Colloid and Interface Science, Vol. 107, No. 1, September 1985