Decoupling of the NernstâPlanck and Poisson Equations. Application

14208

J. Phys. Chem. B 2007, 111, 14208-14222

Decoupling of the Nernst-Planck and Poisson Equations. Application to a Membrane System at Overlimiting Currents Mahamet A.-Kh. Urtenov,† Evgeniya V. Kirillova,‡ Natalia M. Seidova,† and Victor V. Nikonenko*,§ Department of Applied Mathematics, Kuban State UniVersity, 149 StaVropolskaya st., Krasnodar 350040, Russia, Fachhochschule Wiesbaden, UniVersity of Applied Sciences, Kurt-Schumaher-Ring 18, D-65197 Wiesbaden, Germany, and Physical Chemistry Department, Kuban State UniVersity, 149 StaVropolskaya st., Krasnodar 350040, Russia ReceiVed: April 23, 2007; In Final Form: October 9, 2007

This paper deals with one-dimensional stationary Nernst-Planck and Poisson (NPP) equations describing ion electrodiffusion in multicomponent solution/electrode or ion-conductive membrane systems. A general method for resolving ordinary and singularly perturbed problems with these equations is developed. This method is based on the decoupling of NPP equations that results in deduction of an equation containing only the terms with different powers of the electrical field and its derivatives. Then, the solution of this equation, analytical in several cases or numerical, is substituted into the Nernst-Planck equations for calculating the concentration profile for each ion present in the system. Different ionic species are grouped in valency classes that allows one to reduce the dimension of the original set of equations and leads to a relatively easy treatment of multi-ion systems. When applying the method developed, the main attention is paid to ion transfer at limiting and overlimiting currents, where a significant deviation from local electroneutrality occurs. The boundary conditions and different approximations are examined: the local electroneutrality (LEN) assumption and the original assumption of quasi-uniform distribution of the space charge density (QCD). The relations between the ion fluxes at limiting and overlimiting currents are discussed. In particular, attention is paid to the “exaltation” of counterion transfer toward an ion-exchange membrane by co-ion flux leaking through the membrane or generated at the membrane/solution interface. The structure of the multi-ion concentration field in a depleted diffusion boundary layer (DBL) near an ion-exchange membrane at overlimiting currents is analyzed. The presence of salt ions and hydrogen and hydroxyl ions generated in the course of the water “splitting” reaction is considered. The thickness of the DBL and its different zones, as functions of applied current density, are found by fitting experimental current-voltage curves.

1. Introduction The set of the NPP electrodiffusion equations is very important in electrochemistry, being the theoretical background of the mass-transfer description.1-6 The most commonly implicated situation, when the LEN assumption is used instead of the Poisson equation, may be considered as a partial case. From the mathematical point of view, this substitution is possible because the Poisson equation in the dimensionless form contains a parameter (being a factor at the higher derivative) which is generally small. When this parameter is set to be zero, the LEN condition may be applied everywhere in an electrochemical system except at a thin interfacial region (electric double layer) between two different phases, such as the region between an electrolyte solution and an electrode or a charged membrane or interfacial regions in semiconductors. Generally, when applying the LEN condition, the interfacial region is considered as a mathematical plane where the concentrations, electrical potential, and other variables have a stepwise change. Some assumptions are accepted to describe this jump; usually, the local thermo* To whom correspondence should be addressed. E-mail: [email protected]. Phone/Fax: +7 861 219 95 73. † Department of Applied Mathematics, Kuban State University. ‡ University of Applied Sciences. § Physical Chemistry Department, Kuban State University.

dynamic equilibrium is assumed to be valid at the interface, while other conditions are possible.4 However, sometimes it is important to use the Poisson equation essentially. It is the case where a detailed study of the ion transport in the interfacial regions containing a space electric charge is of interest, for example, biological membranes,2,4 the interface between two semiconductors,3 or ion-exchange membranes (bipolar membranes).7,8 Another case is the modeling of intensive concentration polarization in electrode or ionexchange membrane systems where the current density is higher than the classical “limiting” current density. The latter notion is used to define the saturation of the charge transfer by electrodiffusion through a diffusion boundary layer (DBL) not complicated by secondary effects. Following the theories employing the LEN assumption, the interfacial electrolyte concentration in depleting solution tends to zero, and the potential drop through the DBL tends to infinity when the current density tends to its limiting value. The simplest interpretation of the overlimiting current transfer presented in earlier papers9 is the reduction of the effective thickness of the DBL by a growing space charge region. Very important applications of the NPP equations are provided by their coupling with the Navier-Stokes equations. Effectively, the action of the applied electric field on the space

10.1021/jp073103d CCC: $37.00 © 2007 American Chemical Society Published on Web 12/04/2007

Decoupling of the Nernst-Planck and Poisson Equations charge described by the former equations gives rise to an electric body force entering the latter. Thus, mathematically, current passage becomes coupled with liquid flowing. In electromembrane systems, this coupling, known as electroconvection and occurring rather as electro-osmotic slip,10-13 produces the solution mixing near the membrane surface that enhances the mass transfer. The electroconvection phenomenon is also the basis of electro-osmotic micropump functioning14,15 and electrophoretic separation;16,17 the electroconvection coupling is an important element of space charge models describing the ion and water transport in nanometer pores of charged membranes.18,19 Microelectrochemical devices with a thin liquid, solid, or gel electrolyte layer, such as microbatteries,20 present a special case where the Debye screening length is comparable to the film thickness.21 There is a number of papers devoted to the solution of NPP equations. Levich22,23 was the first who formulated the problem and showed that when the current density is less than the limiting current density (Ilim), the DBL near an electrode may be divided into two regions, the electroneutral region in the “outer” part of the DBL and the electrical double layer (EDL) in its “inner” part. The EDL is disturbed (polarized) under the action of the electric current, but it remains at quasi-equilibrium. Grafov and Chernenko,24 Newman,25 MacGillivray,26 and Dukhin et al.27 confirmed this result by means of different mathematical techniques. Smyrl and Newman28 have considered the case of limiting current and found that the region of appreciable deviation from electroneutrality is thicker than the EDL at currents below Ilim and that the potential drop through the DBL including the EDL does not tend to infinity when the current density tends to its classical limiting value. In 1979, Rubinstein and Shtilman9 first found that the NPP equations have a physically meaningful solution, for appropriate boundary conditions, at overlimiting currents in a membrane system. They have shown that at overlimiting currents, the salt counterion concentration reaches its minimum value in a point within the space charge region (SCR) near the interface where the co-ion concentration is negligibly small. With increasing current density, the minimum counterion concentration does not tend directly to zero, as was described before, but decreases rather slightly. Hence, the potential drop through the DBL is still finite while elevated. The local electroneutrality holds within an outer region of the DBL reducing with the current growth while the SCR thickness grows. The increasing concentration gradient in the electroneutral region allows the electrodiffusion to rise over Ilim. Later, a number of numerical solutions of the boundary value problems involving the NPP equations were obtained by Rubinstein et al.,5,12,13,29,30 Mafe´ et al.,31 Manzanares et al.,32 Ben et al.,17 Bazant et al.,21,33and Zabolotsky et al.34,35 Analytical methods for solution of the NPP equation system were developed by Listovnichy36 and Nikonenko et al.,37 who have found some details of the SCR structure by using different approximations resulting from the small parameter method. Urtenov et al.38-40 have found a more rigorous analytical solution for the case of 1:1 electrolyte in a membrane system and formulated39 a new QCD (quasi-uniform distribution of the space charge density) condition that generalized the LEN assumption and allowed one to treat NPP problems by taking into account the deviation from the electroneutrality. Recently, Chu and Bazant33,41 have made an accurate asymptotic analysis for the binary electrolyte case at and above the classical limiting current that allowed them to describe the structure of the depleted boundary layer and to derive a current-voltage relation, which is in good agreement with the numerical solution.

J. Phys. Chem. B, Vol. 111, No. 51, 2007 14209 Zaltzman and Rubinstein30 have developed a unified asymptotic picture of the electric double-layer undercurrent covering all current regimes. A universal electro-osmotic slip formula was deduced from the obtained asymptotics.30 However, all of the reviewed papers deal with particular solutions well-adapted for several cases. Hence, a more universal approach to treat the NPP equations is needed, an approach that could permit one to obtain a rather general algorithm for solving the set of equations in the cases of multi-ion systems, non-steady states, complicated 2D geometries, and so on. Such an approach may be based on the decoupling of the Nernst-Planck and Poisson equations. This envisages a transformation of the initial set of equations in a way that leads to a new system where one of the equations contains only one unknown variable (the electric field) and the others allow the concentrations to be expressed as functions of the electric field. Some aspects of this approach were developed earlier by Urtenov et al.38,39,42,43 Note that Manzanares and Kontturi have given a brief description of the decoupling algorithm in their comprehensive review44 and applied this method for treating several particular multi-ion electrode systems under the LEN assumption. Among others, they suggested to use this method as an effective approximation for evaluating concentration profiles. In complicated cases at underlimiting currents, it is possible to use the Goldman constant field assumption (when the ionic strength is roughly uniform) instead of resolving the differential equation for electric field.44 The purpose of this paper is to develop the decoupling method for the case of a one-dimensional steady-state multi-ion system, which can be applied at under- and overlimiting current densities. We will try to propose a sufficiently general algorithm for obtaining numerical solutions, for asymptotic analysis, and for finding different approximations. The developed method will be used in describing the electrochemical behavior of a multiion membrane system containing two salt (Na+ and Cl-) and two “water” (H+ and OH-) ions, the latter being produced in the course of water splitting at overlimiting currents. Experimental current-voltage curves will be treated and the structure of the depleted DBL analyzed, taking into account the deviation from the LEN assumption. 2. Problem Formulation 2.1. System under Study. As was mentioned above, the set of NPP equations is applicable to different systems having an interface between two phases. In this paper, we will consider a membrane system shown in Figure 1, while sometimes, electrode systems will be discussed also. An electrolyte solution flows between two plane ion-exchange membranes, forming an electrodialysis cell. A direct current passes along the x axis in the direction perpendicular to the membranes. Due to the difference of ion transport numbers in the solution and membranes, the electrolyte concentration decreases near the membranes in one compartment and increases in the neighboring. This current-induced variation in concentration occurring near the interface in a DBL is usually named in membrane science as concentration polarization, while often, this term is regarded in a larger sense as denoting a complex of effects related to the formation of concentration gradients in electrolyte solution adjacent to a permselective (charge selective) solid/ liquid interface upon the passage of an electric current.12,13 Typical current-voltage curves for electrodialysis of dilute solutions are presented in Figure 2. The data are taken from the paper by Zabolotsky et al.45 The potential drop was measured with the help of the Luggin capillaries, the position of their tips being approximately shown in Figure 1. The other characteristics of the system are described in section 10.

14210 J. Phys. Chem. B, Vol. 111, No. 51, 2007

Urtenov et al.

Figure 2. Current-voltage curves for a desalting compartment formed of a cation-exchange membrane MK-40 and an anion-exchange membrane MA-40 with 3 × 3 cm2 of active membrane area and 0.5 mm of intermembrane distance; the positions of the measuring electrodes are shown in Figure 1a. A 0.002 M NaCl solution is flowing between the membranes with a velocity of 3.2 cm s-1. The total and partial currents of the Na+ and H+ ions through the MK-40 membrane are shown. The data are taken from ref 45.

Figure 1. (a) Scheme of an electrodialysis cell with an anion-exchange (AEM) and a cation-exchange (CEM) membrane; DC and CC are desalting and concentrating compartments, respectively. The points 1 and 1′ show where the tips of Luggin’s capillaries are placed. (b) Schematic view of salt counterion (solid line) and co-ion (dashed line) concentration profiles in the DBL of thickness δ near a CEM at different current densities. The DBL is divided in an electroneutral (0 e x e δLEN) and a space charge (δLEN < x e 1) region. The intersection of the straight line extending over the linear concentration profile in the electroneutral region with the x axis gives the effective thickness of the DBL (δ′); λ is the thickness of the quasi-equilibrium boundary layer.

With a growing potential drop, the surface electrolyte concentration in a depleted solution decreases, and the current density of the salt counterion (Na+ ions for the cation-exchange membrane (CEM), in the considered case) is close to being saturated. However, some effects coupled with the passage of intensive currents prevent the saturation.10-13,46-48 The most important one is the electroconvection mentioned above. Besides, gravitational convection46-49 or thermoconvection50 could be involved in several cases, as well as water “splitting”.46,47,51-53 At sufficiently high current densities, the surface concentration of salt ions becomes so small that the H+ (OH-) ions, always presented in solution, start to compete with them in charge transfer, thus giving rise to water splitting at membrane/depleted solution interfaces. The rate of water splitting and, hence, the current density of H+ ions through a CEM (or that of OH- ions through an anion-exchange membrane) depends on the potential difference applied and on the membrane properties. In particular, this rate may be sufficiently high only in the case where functional sites of the membrane are involved in catalytic-type proton-transfer reactions.51-53 2.2. Governing Equations. In this paper, steady-state onedimensional ion transport in a DBL adjacent to an ion-exchange membrane (or an electrode) is considered under the conditions of absent convection. Below, the superscript (u) means that the respective quantity is dimensional (having a unit). We assume that the concentration is uniformly distributed and the local

electroneutrality condition is strictly valid in the center region of channel until x(u) ) 0 (Figure 1). The conditional boundary between the solution and the membrane under study is located at x(u) ) δ(u). Consequently, δ(u) is the thickness of the DBL. The electrodiffusion system of equations for the concentrations ci(u)(x) and the electrical field E(u)(x) within the region between x(u) ) 0 and δ(u) can be written in the form

j(u) i ) -Di

(

dc(u) i

dx(u)

-zic(u) i

dE(u) dx(u)

)

F

F (u) E RT

)

i ) 1, ..., N

(1)

n

∑ zic(u)i (u) i)1

(2)

with the flux densities j (u) i , diffusion coefficients Di, and charge numbers zi of ion i, Faraday’s constant F, gas constant R, absolute temperature T, and dielectric permittivity (u) (which is the product of the dielectric permittivity of vacuum and the relative permittivity of the medium). Note that Kilic, Bazant, and Ajdari54 have recently derived modified NPP equations, which take into account the finite size of ions. This approach is helpful in the case where high concentrations are present (e.g., microbatteries), especially within EDL, where ions become crowded at the surface. However, in membrane systems, the concentrations are rarely higher than 1 mol/L. Hence, the dilute solution approximation, which will be applied in this paper, seems justified. Later on, it is assumed that the ion transfer is stationary in time and no chemical reactions take place in the interval of x(u) ∈ [0, δ(u)]. Under these conditions the fluxes does not depend on x(u). 2.3. Transition to Nondimensional Variables. Let c(u) 10 and D1 be the concentration in the bulk and the diffusion coefficient characteristic for the considered system, where the index “1” will be attributed to the dominant electrolyte counterion. The dimensionless variables are defined as follows

Decoupling of the Nernst-Planck and Poisson Equations

x)

x(u) δ(u)

E)

J. Phys. Chem. B, Vol. 111, No. 51, 2007 14211

(u) c(u) j(u) Di i i δ δ(u) (u) E ci ) (u) ji ) di ) (u) RT D1 c Dc 10

i 10

i ) 1, ..., N I)

I(u)δ(u) D1c(u) 10 F

)

RT(u) (u) 2 c(u) 10 (Fδ )

)

( )

2

L(u) D

2

δ(u)

L(u) D )

x

RT(u)

2 2c(u) 10 (F)

(3)

N

(u) is the current density. where L(u) D is the Debye length and I In electrochemical systems, an important parameter is the limiting current density (I(u) lim). Under this current density, the boundary electrolyte concentration (c1s in Figure 1) becomes small compared to the bulk concentration, and the rate of the potential growth with the current becomes high. In the case where only one kind of counterion is present (species 1) and the membrane is not permeable to co-ions (or only one reaction of discharge of ion 1 occurs at the electrode), we get for a symmetric binary electrolyte1,5

I(u) lim )

2D1c(u) 10 F δ(u)

(4)

(u) where c(u) 10 ) c1 (0) is the bulk concentration of species 1, which is supposedly known. Hence, as it results from eqs 3 and 4, I is equal to the ratio 2I(u)/I(u) lim. By using nondimensional quantities introduced above, eqs 1 and 2 can be rewritten as

dci(x) ) zici(x)E(x) - ji dx

dE dx

2.4. Peculiarities of the Equation System. One important characteristic of the NPP equation system, besides the nonlinearity, is the existence of a small parameter > 0 at the derivative dE/dx. The parameter , which is the squared ratio (u) of the Debye length L(u) D to the DBL thickness δ , usually has a very small value (10-10-10-4). Because of this small parameter presented in the Poisson equation, the boundary value problems for electrodiffusion eqs 5 and 6 are related to singularly perturbed problems.55,56 If is set to be 0, the Poisson eq 6 transforms into the LEN assumption

(5)

F(x) )

zici(x) ) 0 ∑ i)1

(9)

where F is the dimensionless electric space density. Effective analytical methods for resolving eqs 5 and 9 are sufficiently well-known.6,44,57 This problem is sometimes called “degenerated problem”sin relation to a more general problem involving eqs 5 and 6 with * 0. When considering a DBL, the solution of eqs 5 and 9 is restricted by the maximum possible value of the current density that is equal to I(u) lim (eq 4); under this current density, the boundary electrolyte concentration becomes zero. At I(u) > I(u) lim, this concentration becomes negative, which has no physical meaning. For the solution of singularly perturbed equations such as eqs 5 and 6, several general asymptotic and numerical methods have been developed,55-57 not counting solutions,24-40 which are rather particular. However, these general methods are not applicable at overlimiting currents because the solution of the degenerated problem with ) 0 does not exist at I(u) > I(u) lim. Therefore, another general method for the solution of boundary value problems for NPP equations which does not use the degenerated problem is needed. 3. Reduction of the Number of Considered Species. Grouping of Ions into Valency Classes

N

)

zici(x) ∑ i)1

(6)

with x ∈ [0,1] and i ) 1, ..., N. The condition of current passage in dimensionless variables is written as N

diziji ) I ∑ i)1

(7) n1

Below, it will be assumed that the concentrations ci and the electric field E are unknown, whereas the quantities ji and are given as parameters. Equation 7 presents a relation between the fluxes ji (i ) 1, ..., N) and allows one of them, including I, to be found if the others are known. The concentrations in the bulk solution are assumed to be known and linked by the condition of electroneutrality

ci(x ) 0) ) ci0 i ) 1,..., N

The grouping of ionic species was first performed by Schlo¨gl58 in the case of the electroneutrality condition ( ) 0). All of the ions are divided into the groups (classes), and in each of them, only the ions with the same charge are presented. This results in a reduction of the dimension of the equation system. Following Schlo¨gl,58 let us introduce into consideration the “group concentration” Ci and the “group flux” Ji, being the sums of corresponding quantities inside a given group i

(8a)

C1 )

nn

c1j, ...,Cn ) ∑ cnj ∑ j)1 j)1 n1

J1 )

(10)

nn

j1j, ...,Jn ) ∑ jnj ∑ j)1 j)1

(11)

where ni is the number of ions in group i, i ) 1, 2, ..., n, n being the number of groups. By summation of eqs 5 and 6 with the ions of equal charge, we obtain, respectively

N

zici(x ) 0) ) 0 ∑ i)1

(8b)

The particularities of the boundary conditions, especially at the right-hand edge of the DBL (at x ) 1), will be discussed below in section 5.

dCi ) ziCiE - Ji i ) 1, ..., n dx

dE dx

(12)

n

)

ziCi ∑ i)1

(13)


Urtenov et al.

Thus, the dimension of the original set of eqs 5 and 6 is reduced; all of the ion charges now are different, but the form of the equations is still the same. If the number of different ion groups, n, for an electrolyte is 2, we will name this electrolyte “pseudo binary”; n ) 3 will correspond to a (pseudo) ternary electrolyte, and so forth. 4. Decoupling of the Nernst-Planck and Poisson Equations As was mentioned above, eqs 5 and 6 (and, evidently, eqs 12 and 13) belong to singularly perturbed nonlinear ordinary differential equations, which are generally difficult to solve. From the structure of the equations, it follows that if the electrical field is determined, the concentrations Ci (i ) 1, ..., n) can be found easily by the solution of the Nernst-Planck eq 12, which is linear with respect to Ci. Note that this approach was proposed, according to Manzanares and Kontturi,44 by Kramers.59 In this paper, we will derive an equation for E which is independent of the concentrations Ci. Therefore, all unknown quantities E and C1, ..., Cn can be determined separately from independent equations. Let

n

a() )

n

1

Ci0 - E20 ) ∑ Ci0 + O() ∑ 2 i)1 i)1

Here, Ci0 ) Ci(x ) 0). The electric field at x ) 0, E0, is on the order of 1, as can be seen from eq 5. Note that the integration uses the fact that under steady-state transfer and in the absence of chemical reactions, the fluxes are constant. Multiplication of eq 12 by zki (k ) 0, 1, ..., n - 1) and subsequent summation from i ) 1 to n yields a recurrence formula

d S ) Sn,k+1E - Gn,k k ) 0, 1, ..., n - 1 dx n,k

1 Sn,0 ) E2 - (Gn,0x - a) 2 Sn,1 )

zki Ci(x, ) ∑ i)1

n

Gn,k )

zki ji ∑ i)1

k ) 0, 1, ..., n (14) Sn,2 )

Note that n

Sn,0 )

Ci ∑ i)1

n

Sn,1 )

ziCi ∑ i)1

Sn,3 )

n

Sn,2 )

z2i Ci ∑ i)1

are, in dimensionless notations, the total concentration, the charge density, and the double ionic strength, respectively. The use of linear combinations of concentrations and fluxes was apparently introduced by Schlo¨gl.58 The advantages of using the ionic strength as a variable in mathematical modeling of ionic transfer in multicomponent solutions were admitted by other authors.60 Sn,k and Gn,k can be considered as “supraconcentrations” and “suprafluxes”, respectively. Next, we introduce quantities qn,k, which are the sums of possible products of k terms over z1, ..., zn, for example

q4,1 ) z1 + z2 + z3 + z4 q3,2 ) z1z2 + z1z3 + z2z3 q3,3 ) z1z2z3 (15) From the definition of Sn,k (k ) 0, 1, ..., n), it follows that Sn,n can be expressed by Sn,0, ..., Sn,n-1

S2,2 ) q2,1S2,1 - q2,2S2,0

(16a)

S3,3 ) q3,1S3,2 - q3,2S3,1 + q3,3S3,0

(16b)

S4,4 ) q4,1S4,3 - q4,2S4,2 + q4,3S4,1 - q4,4S4,0

(16c)

We find the sum of eq 12 for all i and insert eq 13 into the result. The integration of the obtained equation from x ) 0 to an arbitrary point x leads to the first integral n

Sn,0 )

1

Ci ) E2 - (Gn,0x - a) ∑ 2 i)1

where the integration constant a is expressed as38

(17)

(19)

Thus, the system of n Nernst-Planck equations is replaced by n equations under the form of eq 19. The functions Sn,0, Sn,1, ..., Sn,n-1 can be determined now successively and expressed by the function E and its derivatives. Really, Sn,0 is already found (eq 17), and its substitution into eq 19 gives Sn,1, and so on

Sn,k(x,) ) n

(18)

dE dx

(

(20b)

1 d2 E + Gn,1 E dx2

[ (

(20a)

)

(20c)

) ]

Gn,2 1 d d2E Gn,1 + + E dx E dx2 E E

(20d)

Equations 15-20 are presented for (pseudo) binary (n ) 2), ternary (n ) 3), and quaternary (n ) 4) electrolytes. Usually, these systems are sufficient for the treatment of the majority of practical cases. Note that eq 20a for a binary electrolyte was first derived by Grafov and Chernenko.24 To deduce an equation containing only the electric field E and its derivatives, eq 19 is written for a given n and k ) n 1. Then, Sn,n-1 on the left side of eq 19 is substituted by a corresponding eq 20. Sn,n present on the right side of eq 19 is expressed in terms of Sn,n-1, Sn,n-2, and so forth by using eq 16; the latter are then expressed in terms of E by eq 20. The following equations for E are obtained, respectively, for n ) 2 and 3

d 2E ) dx2 q2,1E

[

dE - q E3 + q2,2(G2,0x - a)E - G2,1 (21a) dx 2 2,2

]

[

]

G3,1 dE d2E dE d 3E ) + + q3,1E - q3,2E2 + 3 2 E dx E dx dx dx q E4 + q3,3(a - G3,0x)E2 + G3,1q3,1E - G3,2E (21b) 2 3,3

A similar equation is deduced for the case n ) 4.43 In the case of a binary electrolyte, eq 21a, which is the master equation for the electric field, was first derived by Smyrl and Newman,28 while Grafov and Chernenko24 obtained a single differential equation for the co-ion concentration. Urtenov et al.,38,39,61 Rubinstein and Zaltzman,12,29,30 and Chu and Bazant33,41 have made the detailed analysis of this equation at different current densities.



The list of equations for eq 21 for n greater than 3 might be continued, but their form becomes rather complicated. Equations 21 can be solved numerically by the finite difference method together with Newton’s and Thomson’s methods,62 while other methods are also used.33 They are also suitable for asymptotic analysis.

E

|

x)1

)

x [∑

x

2

1

5. Boundary Conditions

q2,2(G2,0x - a)E - G2,1 ) 0

(22)

The solution to eq 22 is n

E(x) )

(

G2,1

q2,2 G2,0x - a

Ci(1) + Gn,0 - a )

i)1

x

According to the theory of singularly perturbed problems,55-57 eq 21 with an unknown electric field and known fluxes demands boundary conditions at both limits of the interval x ∈ [0,1]. Usually, it is assumed that the bulk concentrations ci0 are known (eq 8a) and satisfy the electroneutrality condition (eq 8b). Then, the value of the electric field at x ) 0 can be found. Consider, for example, the case of the (pseudo) binary electrolyte. As was mentioned above, eq 9 and, hence, eq 8a are obtained from the Poisson eq 6 when setting ) 0. Equation 21a under the LEN assumption (with ) 0) becomes

zi Ji ∑ i)1

) (∑ )

n

|z1z2|

i)1

n

Ci0 - x

∑ i)1

)

(23)

Ji

For x ) 0, we have E(0) ) -G2,1/aq2,2. It can be seen that E(0) is on the order of 1, which justifies eq 18 for the integration constant a. Substituting eq 18 (without omitting the small term 0.5E20) and taking into account the definitions (eq 14) for Gn,k gives the last expression in eq 23. The construction of the boundary conditions at x ) 1 is more difficult because of the undistinguished boundary between the solution and solid phases. This problem arises from the complicated structure of the interface and the vagueness in the question of under which conditions the NPP equations are applicable;5,6,54,63 for example, what is the distance from a boundary atom (molecule) belonging to the membrane starting from which the NPP equations (in the form of eq 5) can be used in the solution? In this paper, we leave this problem. However, for the sake of concreteness, we assume that x ) 1 relates to the inner edge of the diffuse part of the electrical double layer, that is, the outer edge of the compact Helmholtz (or Stern in the case of adsorption) layer. If the electroneutrality condition is used and the fluxes Ji are known, no boundary conditions at x ) 1 (which corresponds, in this case, to the outer edge of EDL) are needed; the concentrations and the electric field in the solution at any x ∈ [0,1] are found from the Nernst-Planck equations if their values at x ) 0 are given. To find the relation between the jumps of concentrations and the potential across the interface, the Donnan local equilibrium equation may be applied to the interface considered as a plane.2 In the general case where the Poisson equation is applied and the equations of eq 21 are to be solved, the electric field should be given at x ) 1. As follows from the first integral, eq 17, E(1) should be in accordance with the concentrations at x)1

]

n

1

x[

n

2 Gn,0 +

[Ci(1) - Ci(0)] ∑ i)1

]

(24)

The term 0.5E2(0) is omitted in the last expression on the righthand side of eq 24. Ci(1) and Gn,0 can be considered as parameters of a given boundary value problem. Their values may be inserted as known or found from other conditions. For example, to find Ci(1), note that, as follows from the Poisson eq 6, E(1) is equal to the dimensionless electric space charge in the solution

E(1) )

∫01 ∑ ziCi(x)dx ) ∫01 F(x)dx i

Then, equations similar to eq 24 can be deduced in the phase neighboring the solution, for example, an ion-exchange membrane. By equating the overall space charges in both phases (the condition of global electroneutrality) and using an appropriate condition for linking the boundary concentrations, an equation for determining Ci(1) may be found. Instead of E|x)1, the derivative (dE/dx)|x)1 can be given. Its value should be in accordance with Ci(1) also

dE dx

|

F(1) ) x)1

)

1

∑i ziCi(1)

(25)

It will be shown (section 9) that the use of a new assumption named the assumption of quasi-uniform charge density distribution (QCD) does not demand any condition on concentration or electric field at x ) 1. The boundary conditions in a system consisting of two parallel electrodes and a thin-film electrolyte between them were analyzed by Bazant et al.,21,33,54 by Zaltzman and Rubinstein,30 and by other authors cited therein. They used an integral constraint condition, reflecting the fact that the total number of anions within the film (which do not participate in electrochemical reactions and are not subject to specific adsorption) is constant. The authors21,33,54 considered the compact part of EDL and link the rate of the Faradaic electron-transfer reaction with the potential drop over the Stern layer by the Butler-Volmer equation, following the Frumkin approach. Then, specifying the current (which is transported exclusively by cations), they could resolve the NPP equation system. 6. Numerical Solution The presence of small parameter in eq 21 significantly complicates the integration. It is well-known9-13 that, in this case, the electric field E and concentrations Ci vary very rapidly in a narrow boundary layer, the thickness of which is close to the Debye length L(u) D . The dimensionless analogue of the Debye length is on the order of x. One of the possible solutions for this problem is the use of nonregular grid condensing near the interface.32 Another way is the use of new variables38,61 -2/3 -1/3 ξ ) (Gn,0x - a)Gn,0

(26a)

-1/3 1/3 η(ξ,) ) E(x,)Gn,0

(27b)


Urtenov et al.

Figure 3. The distributions of dimensionless electric fields E(x) (a), total cation and anion concentrations (b), counterion Na+ (cNa) and co-ion Cl(cCl) and OH- (cOH) concentrations (c), as well as the space charge density F(x) (d) in the DBL near a cation-exchange membrane found by numerical solution of NPP equations for the following dimensionless parameters defined by eq 3: j1 ) jNa ) 2.979, j2 ) jCl ) 0, j3 ) jH ) 0, j4 ) jOH ) 0.577, and ) 4.56 × 10-7. The parameters correspond to a membrane system described in the text and in Figure 2 at I(u) ) 9.52 mA (u) cm-2, ∆φ(u) ) 14.2 µm. The calculations are made for three different values of c1(1) as a boundary condition at x ) 1, 0.1, 1, tot ) 7.7 V, and δ and 10 (the values are indicated in the legends), and under the QCD assumption not requiring c1(1). The vertical straight line on the left of δ′ shows the value of δLEN ) δ′ - 1/3, and the vertical straight line on the right of δ′ shows the inner edge of the transition zone calculated as δ′ + 1/3.

Equation 21a is transformed then into eq 27

G2,1 dη 1 d2η ) q2,1η - q2,2η3 + q2,2ξη 2 dξ 2 G2,0 dξ

(27)

Similar equations are obtained for n ) 3 and 4.61 Equation 27 does not contain small parameters; eq 26a stretches the argument coordinate axis, and eq 26b reduces the function axis. As a result, the derivative dη/dξ does not attain high values typical for dE/dx, which allows for numerical integration of eq 27 with a regular grid for any values of . The distribution of the electric field E(x) resulting from the function η(ξ) as well as the group and individual ion Ci(x) and charge density F(x) are shown in Figure 3a, b, c, and d, respectively. The depleted DBL near a cation-exchange membrane is considered. Four ionic species are present, two cations (Na+ and H+) and two anions (Cl- and OH-); the concentration of the H+ ions is linked with that of the OH- ions by the relation CHCOH ) KW, where KW is the ion product constant of water. There is 0.002 M NaCl in the bulk (x ) 0). The OH- ions generated in the course of water splitting at the membrane interface (x ) 1) occurring at I(u) > I(u) lim are present only inside of the diffusion layer. The electrolyte is treated as pseudo binary. In calculations, we used the parameters relating to the membrane system, the current-voltage curves of which are shown in Figure 2. The total current density, I(u), was taken to be equal

to 95.2 A m-2, and the partial current densities of the Na+ and H+ ions through the MK-40 membrane were 54 and 41.2 A m-2, respectively (as measured in ref 45), the latter being equal to the OH- ion current density in the DBL; the currents of Cland H+ in the DBL were assumed to be zero. The DBL thickness (δ(u) ) 14.2 µm) was fitted in such a way that the potential drop across it, calculated by integration of E(x), was close to that estimated from the experiment (≈4 V) (see also section 10). Three values of the dimensionless counterion concentration at x ) 1 (cNa(x ) 1)) were taken, 0.1, 1, and 10. Analysis of the curves presented in Figure 3a, b, c, and d is given in the following sections. 7. LEN Assumption. Limiting Current 7.1. The Condition of the Limiting State for a Pseudo Binary Electrolyte Solution. As was mentioned above, the solution of the NPP equations in the case of the LEN assumption is obtained when the small parameter is set to be zero, ) 0.1,5,33 In this case, the solution of eq 21a obtained from the NPP equations for a (pseudo) binary electrolyte leads to eq 23. The substitution of eq 23 into eq 12 and the subsequent integration leads to linear concentration profiles of the total group concentrations of cations (C1) and that of anions (C2). When a CEM is considered and I > 0, the parameters G2,0 and G2,1 are positive, q2,2 ) z1z2 < 0. Then, eq 23 has a physical



meaning if (G2,0x - a) < 0. To satisfy this condition for all x ∈ [0,1], it is necessary that

G2,0 < a a ≈

∑i Ci0

(eq 18)

Thus, in the case where the LEN assumption is applied, the supraflux G2,0 cannot possess values equal to or higher than a. If

G2,0 ) G2,0 lim ) a

(28a)

E(1) ) ∞; this state corresponds to the limiting current density. In the limiting state, the concentrations of all ions at x ) 1 are zero. In dimension quantities, the condition of the limiting state G2,0 ) a is written as n1+n2 n1+n2

∑ i)1

(u) j ilim

)

Di

∑ i)1

(u) ci0

(28b) δ(u)

if the term O() is neglected in eq 18, expressing a, n1 and n2 being the numbers of cations and anions, respectively, in the pseudo binary electrolyte. 7.2. Limiting Fluxes, Limiting Current Density, and the Exaltation Effect. Consider a particular case where n1 ) n2 ) 1, one cation (i ) 1) and one anion (i ) 2) are present in the solution. If a cation-exchange membrane is considered, these ions are a counterion and a co-ion, respectively. Then, the flux density of the counterion can be expressed from eq 28b as follows

j (u) 1lim )

2D1c(u) 10 δ

(u)

-

D1 (u) j D2 2lim

(29)

The charges z1 and (-z2) are set as 1, for the sake of simplicity. The transfer of the co-ion through the membrane may be characterized by its effective6 (or integral44) transport number (u) T2 ) z2j(u) 2 F/I , which is a measure of the membrane nonideal permselectivity. By using this definition and the fact that the current density is

I(u) )

∑i zij(u)i F

the following well-known expression for the limiting current density may be found from eq 29

I (u) lim

)

Dc(u) 0 F (t2 - T2)δ(u)

(30)

Note that eq 30 can be also obtained when the charges zi are arbitrary, with D ) (z1 - z2)D1D2/(z1D1 - z2D2), the diffusion (u) coefficient of the electrolyte; c(u) ) z1c(u) 0 10 ) -z2c20 , its -1 concentration (in eq L ); and t2 ) |z2|D2/(z1D1 - z2D2), the co-ion transport number in solution. If the flux of co-ions, j(u) 2 , is zero, eq 30 transforms into eq 4. Consider also a case where two kinds of co-ions are present in the depleted diffusion layer. The first co-ion (A-) belongs to a salt electrolyte MA. The second one may appear there due to an electrochemical reaction. In membrane systems, water splitting usually occurs at the membrane/solution interface under

currents close to or higher than the limiting one.45-52 If the membrane is a cation-exchange one, H+ ions pass through the membrane, whereas OH- ions move toward the solution bulk. In an electrode/solution system, a similar situation is possible when a second electrode process produces co-ions, for example, when there is electroreduction of oxygen, injected OH- ions into the solution, at a metal electrode in parallel with the reduction of a cation M+.64,65 Let j2(u) be the flux density of coions (OH- ions) generated in an electrochemical reaction. Then, the limiting salt counterion flux can be calculated by eq 29, similar as that for a 1:1 electrolyte, with the only difference being that index “2” should relate now to OH- ions, the salt co-ion flux being zero. Thus, eq 29 shows that regardless of the mechanism for coion flux occurring in the DBL, the flux density of the counterion depends on it. In the limiting state, according to eq 29, the salt counterion flux can be presented as a sum of three components, a diffusion contribution, a migration one under an electric field assumed to not be disturbed by co-ion transfer (i.e., when setting j2 ) 0), and an additional component known, in the case of the parallel electrochemical interfacial reaction, as the “exaltation current”.64-67 Note that in the case of a symmetrical salt electrolyte, the first and the second components are equal to (u) (thus giving the first term on the right-hand side of D1c(u) 10 /δ (u) eq 29), the third component j (u) 1exalt ) -(D1/D2)j 2lim being (u) positive as j2 is negative in the considered case. The reason for exaltation is an increase in the electric field produced by co-ions penetrating through or generating at the membrane (electrode)/solution interface. This increase is easily seen from eq 23; in the limiting state, G2,0 ) a, whereas G2,1 (which is the dimensionless current density) and, hence, E(x) grow with (u) an increase in | j (u) 2lim|. In other words, the growth in j 1lim with (u) increasing | j 2lim| may be explained by the fact that the co-ions appearing near the interface attract the counterions from the solution bulk. In membrane science, eq 29 as applied to describe the exaltation effect is known under the name of the “Kharkats equation”, in honor of this scientist who first studied the phenomenon.66 Equation 29 shows clearly that the exaltation of the counterion flux occurs by the co-ions independent of the nature of the latter and the reason why they are appearing at the membrane (electrode)/solution interface; it can be a parallel electrochemical reaction generating, for example, OH- ions as considered above or the leakage of salt co-ions through the membrane. It is interesting to note that the counterion flux exaltation by salt co-ions was never mentioned in the literature, as far as we know, while eq 30 or similar equations68 linking the current density and the co-ion flux are known. Traditionally,69 when deducing eq 30, the diffusion and migration contributions to the ionic fluxes are considered, and this equation is obtained by equaling the fluxes in the DBL and the membrane. The limiting current density is then reached when the diffusion contribution attains its maximum. The fact that the migration and, hence, the overall flux of counterions increases with increasing co-ion leakage through the membrane seems unnoticed. 7.3. Competition between Two Counterions. In the case where two different competitive counterions (e.g., Na+ and K+) and a co-ion (Cl-) are present in the solution, the limiting flux density of one of them is (u) j 1klim )

(1 - z1/z2)D1kc(u) 1k0 δ

(u)

+

z1D1kc(u) 1k0 z2D2c(u) 20

j (u) 2lim k ) 1, 2

(31)


Urtenov et al.

where index “1k” is reserved for the counterion under consideration entering ionic group 1, and index “2” is for the co-ion group. In this case, the concentration profiles of the counterions are linear. Different counterions (Na+, K+) behave as explicitly independent each of others; however, both “feel” the overall flux of co-ions through the membrane. If the bulk concentration of a counterion with “12” index decreases, the co-ion concentration c(u) 20 becomes lower also, and this results in an increase of the co-ion influence on the “11” counterion, as it follows from the exaltation term of eq 31. 7.4. Pseudo Ternary Electrolyte. For a (pseudo) ternary electrolyte, similar conclusions can be obtained. In this case, by omitting in eq 21b all terms with , we come to an Abelian equation of the first kind

(

)

G3,2 2 dE q3,3 ) (G3,0x - a)E3 - q3,1 E dx G3,1 G3,1

(32)

E(x) can be found by resolving a Cauchy problem62 for eq 32 with the boundary condition n

/∑

E(0) ) G3,1

z2i ci0

(33)

i)1

It can be shown that the limiting fluxes of competitive counterions can be calculated by the same eq 31 as that for a binary electrolyte. Hence, the shape of this equation does not depend on the number of different counterion types; the charges can take arbitrary values.

boundary. The solution of eq 21a without this term (the QCD assumption, section 9) is shown in Figure 3a with a dotted line. 8.2. Concentration and Space Charge Density Distribution. Figure 3b shows the profiles of summary cation (cNa + cH) and anion (cCl + cOH) concentrations, and Figure 3c shows the concentrations of three ions presenting in the depleted DBL one counterion, Na+, and two co-ions, Cl- and OH-; the concentration of H+ is small and not shown. The profiles are calculated numerically from the Nernst-Planck eq 5 after substituting the electric field calculated previously. Four zones of the DBL can be distinguished, as well as in the case where only two ions are present.9,37,38 The first one is the zone of local electroneutrality extending from x ) 0 to δLEN, eqs 9 and 23 being verified there. The second one [δ′ + λtr, xs] is the migration37,38 (or extended space charge)30 zone of the SCR, where the concentrations vary only slightly and the main mechanism of ion transport is electrical migration, the contribution of diffusion being negligible. In the third zone, the gradients of the counterion concentration and the electric field are both high. The first and second zones are separated by a transition zone, the thickness of which is on the order of 1/3.28,30,33,37,38,61 In the third zone, the concentration of co-ions is negligible, and the diffusion and electrical forces acting on the counterion are directed to the opposite sides. Thus, their gradients are very high, but the flux is small. As a consequence, the flux term for the counterion (j1) can be neglected in this zone compared to the diffusion and migration terms of the Nernst-Planck eq 5 in the leading-order approximation.33,38,39 Thus, the relation between the functions c1(x) and E(x) is the same as that when the current is zero; it is the Boltzmann equilibrium distribution

8. Overlimiting Current 8.1. Electric Field Distribution. The curves E(x) presented in Figure 3a were obtained by numerical solution of eq 21a at three different values of the electric field at the right boundary E(1): 1485, 2481, and 6756 or 2.68, 4.48, and 12.2 (106 V m-1) when passing to the dimension quantities. These values of E(1) are calculated from eq 24 after setting C1(1) ) 0.1, 1, and 10, respectively. Figure 3a shows also the electric field calculated using eq 23, corresponding to the LEN assumption. The fact that the current density surpasses the limiting one (I ) 5.25 in calculations whereas Ilim ) 2) signifies that, in the considered case of a binary solution, G2,0 > a ≈ 2, effectively, G2,0 ) 2.40. The condition G2,0x > a necessary for the remaining LEN eq 23 to be valid is satisfied for 0 e x < δ′, where δ′ ) a/G2,0. δ′ is the effective thickness of the DBL; it can be obtained by the cross section of the salt concentration profile with the x axis (Figure 3b). When x f δ′, E calculated by LEN eq 23 tends to the infinity (Figure 3a). However, in the region of 0 e x < δLEN, where δLEN is a point lying at a small distance λtr (of the order of (/I)1/3)38,61 on the left of x ) δ′ (Figures 1 and 3a,b), E found from eq 23 coincides with that calculated numerically by resolving the “exact” eq 21a. Hence, the electroneutrality assumption can be used in the interval [0, δLEN] for calculating E. However, at this point, the concentration of the co-ion differs significantly from that of the counterion, and the space charge density is far from zero (Figure 3b and d). It is important to note that the value of the electric field at the right boundary E(1) affects the E(x) curve only in a narrow zone near the interface. Outside of this zone, E(x) does not depend on E(1). The explanation is that the term d2E/dx2 on the left side of eq 21a is significant only in the boundary zone called quasi-equilibrium (see the next subsection). In the other regions, this term can be omitted, and eq 21a becomes of the first order; hence, it does not demand conditions at the right

dc1(x) ≈ z1c1(x)E(x) dx

(34)

Hence, the third zone may be called the quasi-equilibrium space charge zone30,37-39 (it was named the “inner diffuse layer” by Chu and Bazant).33 Its thickness (λ) is on the order of the Debye length determined by the counterion concentration at the outer limit of this zone, namely, by its minimum concentration (c1s) in the DBL (Figure 1b), λ∼ x/c1s.37-39 When the current increases, c1s decreases starting from c10 at I ) 0 and is on the order of 1 at I < Ilim; c1s becomes on the order of x37-39 at I > Ilim. Hence, λ is on the order of 1/2 when I < Ilim and on the order of 1/4 when I > Ilim. It follows from eq 17 that in the electroneutral zone, Sn,0 depends linearly on x, as can formally be taken as equal to 0 in this zone. In the considered case, C1 ) cNa + cH ) cCl + cOH ) 0.5Sn,0; hence, cNa + cH as well as cCl + cOH are linear (Figure 3b). If the salt co-ion flux is zero, the concentration of this ion can be expressed by the following equation of the Boltzmann type obtained by integration of the Nernst-Planck eq 5

| |∫

c2 x ) c02exp - z2

()

(

x

0

||

E(l)dl ) c02exp - z2 ∆φ(x)

)

(

)

(35)

where ∆φ(x) is the potential drop in the interval [0, x]. In the electroneutral zone, ∆φ(x) grows with the rate of a logarithm, as can be seen from eq 23; hence, c2(x) decreases with x (linearly in the case cOH ) 0). In the space charge region, E(x) grows strongly (Figure 3a) but does not tend to infinity; hence, c2(x) decreases with x very rapidly. Consider now the transition zone. Its thickness, as was mentioned above, is on the order of 1/3. This value is normally small in comparison with the entire DBL; as well, the electric



field, E, in this zone is relatively low at overlimiting currents when the maximum of E falls on the migration zone. Hence, when calculating approximately the overall voltage at I > Ilim, the peculiarities of this zone may be ignored.37 However, the transition zone is important for understanding the behavior of the electrochemical system, especially at I ) Ilim. The existence of this zone at I ) Ilim was first discovered and its behavior studied by Smyrl and Newman.28 The transition zone arises at current densities lower than but very close to Ilim (depending on ) when the salt boundary concentration (cs ) c1s ) c2s) becomes on the order of 1/2. Then, with increasing current, c2, at the distances on the order of 1/3, continue to decrease, while c1s decreases only slightly (Figure 1b). At higher currents, the outer edge of the transition zone corresponds to a point where the deviation from the electroneutrality becomes obvious (Figure 3b); its inner edge joins the migration zone, the thickness of which increases with I starting from 0 at I ) Ilim.37 In the transition zone, the rate of decrease in concentration of the counterion with x is lower than that of the co-ion. As a result, the space charge density F(x) increases rapidly in this zone, being on the order of in the electroneutral zone and O(1/2) in the transition one;30,33,38,39 see also the next section. However, when x continues to grow, c1(x) continues to decline, whereas c2(x) becomes negligible. Thus, F(x) has a local maximum within the transition zone (Figure 3d) that is at a distance close to (1 - δ′) from the surface. As the action of the viscosity forces gets weaker with increasing distance from the surface, the possibility to generate fluid flow is higher when the volume force is applied farther from the surface. This fact is significant when explaining the phenomenon of electroconvection occurring as electro-osmotic slip along the surface.11-13 The exact limits of the transition zone may be defined only approximately as the difference in concentrations of cations and anions increases monotonically with increasing x at the beginning of the SCR. The interval [δ′ - 1/3, δ′ + 1/3] is rather a lower estimation of this zone. It is important to observe that the maximum of F lies inside of this zone, at least under the conditions of the calculations made in this study. The fact that the DBL is distinctly divided into electroneutral and space charge regions allows us generalization of eqs 28 and 29 to the case of the overlimiting currents. Effectively, when treating eqs 5 and 6 in the region of 0 e x e δ′ at I > Ilim in the same way as earlier in the overall DBL (0 e x e 1) at I ) Ilim, we find an equation similar to eq 29

j (u) 1 )

2D1c(u) 10 (u)

δ′

-

D1 (u) j D2 2

(36)

where δ′(u) is the effective thickness of the diffusion layer; δ′(u) is equal to the distance from x ) 0 to the point where the straight line following the counterion concentration profile in the electroneutral region cuts the x axis (Figure 1b). (Remember that the superscript (u) is used to show that the quantity is dimensional, δ′ ) δ′(u)/δ(u)). Similarly, it can be shown that eq 31 for the fluxes of competitive counterions with an arbitrary charge is verified at I > Ilim; if δ(u) is substituted by δ′(u), the number of counterions has no importance. 9. QCD Assumption As follows from the numerical solution, the term (d2E/dx2) in eq 21a can be neglected in the region of [0, xs], where xs )1 - λ; λ is the dimensionless thickness of the quasi-equilibrium part of the SCR, and xs relates to the point where the counterion

concentration attains its minimum (Figure 1b and c). Consider the problem of the approximate solution of eq 21 in a more consistent way. When x increases from 0 to xs, E(x) increases from E0 ) E(0) ) O(1) to Es ) E(xs). At overlimiting currents, the minimum value of the salt counterion concentration c1s achieved at x ) xs is small compared to 1, being on the order of x.38,61 The concentration of co-ions at this point is much smaller. Hence, the left-hand term of the first integral, eq 17, is much smaller than 1, being on the order of x. Then, as can be seen from eq 17, the term (/2)E2s should be on the order of 1 because the second term in the right-hand part of this equation, (Gn,0xs - a) is on this order. Thus, we find that Es is on the order of 1/x. As E(x) is a monotonous function (Figure 3a), its derivative in the interval of [0, xs] can be estimated as dE/dx ∼ Es/xs ∼ 1/ x, xs ∼ 1, that is, dE/dx ) O(1/x). Similarly, the following estimations of the terms of eq 21a can be made

( )

d2 E dE 1 ) O(x) E ) O(1) E3 ) O dx dx2 x

Hence, the term (d2E/dx2) in eq 21a has the smallest order of magnitude for x ∈ [0, xs] when parameter takes small values (normally, ∼ 10-10 to 10-6). This term can be neglected by taking into account that qn,k and Gn,k are on the order of 1. Therefore, in the case of a (pseudo) binary electrolyte and x ∈ [0, xs], we obtain

q2,1E

dE ) q E3 - q2,2(G2,0x - a)E + G2,1 dx 2 2,2

(37)

For a symmetric electrolyte with z1 ) -z2 (q2,1 ) 0), the cubic equation

q E3 - q2,2(G2,0x - a)E + G2,1 ) 0 2 2,2

(39)

is the master equation for the electrical field E. Equation 38 is written as

3 E - (Ix - a)E - I ) 0 2

(30)

for an ideally selective membrane (j2 ) 0) and a 1:1 electrolyte; q2,2 ) -1, G2,0 ) G2,1 ) I. The passage from eq 21a to eq 37 may be interpreted as follows. It can be seen from the Poisson eq 6 that

d2E dF ) dx2 dx

(40)

When setting ) 0 and obtaining eq 22 from eq 21a, we assume effectively that F ) ∑ zici ) 0, that is, the LEN condition. When omitting only the term (d2E/dx2) in eq 21a and obtaining eq 37, we assume that

d F(x,) ) 0 dx

(41)

It means that the function F(x,) is no longer negligible, but its derivative (d/dx)F(x,) can be omitted in the region of [0, xs], being on the same order of magnitude as x. The latter follows from eq 6 and the estimation of Es made earlier. Equation 41 can be interpreted also in such a way that the charge density is assumed to be distributed nearly uniformly in


Urtenov et al.

the region of [0, xs]. It is the condition of quasi-uniform charge density distribution (QCD). In the general case (I < Ilim or I g Ilim), the LEN and the QCD conditions have their area of application; the electroneutrality condition can be applied at x ∈ [0, δLEN] and the QCD condition at x ∈ [0, xs]. Note that the LEN condition does not strictly mean that F ) ∑ zici ) 0 in [0, δLEN]; really F is on the order of there, as it follows from the Poisson eq 6 and the fact that E(x e δ1) is on the order of 1. At x ∈ [δLEN, xs], the order of magnitude of F is the same as that of c1s, that is, x; that of (d/dx)F(x,) is also x, as was mentioned above. In the quasi-equilibrium zone of SCR [xs, 1], the order of derivative d/dx F(x,) becomes of the order of -1/4 (as dF ∼1 and dx ∼ λ ∼ x/c1s ∼ 1/4); thus this term cannot be omitted in eq (21a). While the QCD condition is not applicable in the quasiequilibrium zone, it is possible to find the electric field and concentration distribution in this zone. An approximate eq 34 and eq 6 with F ≈ z1c1 can be used for this aim. The potential drop over the quasi-equilibrium part of SCR can be found by the integration of eq 34 from x ) 1-λ to 1 applied to a counterion

∆φsD )

1 1 E(x)dx ≈ ln ∫1-λ z1

c1(1) c1s

(42)

Note that eq 34 can be used not only in the quasi-equilibrium interfacial layer of the solution but in that of the membrane as well. As there are fixed ions with relatively a high concentration (C h X) of about 1 M, the SCR thickness there (λh) is very small, on the order of the Debye length determined by the C h X value, and only slightly depends on the applied current density. Under these conditions, eq 34 can be integrated over the interval of [1 - λ, 1 + λh], and instead of eq 42, the following equation is obtained

∆φID )

cj

1+λh 1s 1 E(x)dx ≈ ln ∫1-λ z1 c1s

(43)

where cj1s is the counterion “surface” concentration in the electroneutral bulk membrane at the interface (I) with the diluting solution (Figure 1). Index D is used because eq 42 (as well as eq 43) is the Donnan equation for the potential drop over an equilibrium interface. The estimation of the order of magnitude of the terms in eq 21b in the case of the ternary electrolyte is similar to that presented above for a pseudo binary electrolyte. In this case, one obtains a master differential equation valid at x ∈ [0, xs], which is on the same first order as that of eq 32, found when applying the LEN assumption.61 A similar situation takes place for electrolytes with a larger number of species.61 Thus, the application of the QCD assumption to one of the equations in eq 21, written for the electric field, permits one to simplify it and to obtain a first-order differential equation or an algebraic equation (in the case of a symmetric pseudo binary electrolyte). The solution of the obtained equation with an appropriate condition at x ) 0 gives the E(x) function everywhere (including the electroneutral, the transition, and the migration zones) except the quasi-equilibrium zone of SCR with the thickness λ on the order of x/c1s near x ) 1. The substitution of the found E(x) into the Nernst-Planck equation written for species i leads to a Cauchy problem, whose solution yields the concentration distribution ci(x). The QCD assumption does not demand any conditions at x ) 1 if the fluxes of all

ions are given. The counterion concentration calculated in this way at x ) 1 and the overlimiting currents is very close to c1s, for example, only 2% lower than c1s, under the conditions of Figure 3c. The potential drop over the quasi-equilibrium zone can be found by the equilibrium Donnan eq 43 without calculation of the E(x) and ci(x) values inside of this zone. Thus, the QCD assumption agrees with the condition of local equilibrium at the interface, which can be applied to an interfacial layer with the thickness on the order of λ + λh. In electrodialysis membrane systems, this layer is normally very thin compared to the DBL thickness. Thus, for calculating the voltage-current curves, the interfacial layer can be considered as a Donnan quasi-equilibrium plane. Then, the QCD assumption allows us to obtain a boundary value problem for under- and overlimiting currents with the same boundary conditions as those when using the LEN assumption. 10. Treatment of Experimental Current-Voltage Curves for a Membrane System. Discussion 10.1. Membrane System. The current-voltage characteristic of an electrodialysis desalination channel studied in ref 45 is examined. The channel is formed of two heterogeneous membranes, a cation-exchange MK-40 and an anion-exchange MA-40 membrane produced by SchekinoAzot, Russia. The active membrane area is 3 cm × 3 cm, and the intermembrane distance 0.5 mm. A 0.002 M NaCl solution is flowing between the membranes, with an average velocity of 3.2 cm s-1. The total and partial current densities of the Na+ and H+ ions through the MK-40 membrane found by using a special technique45 are shown in Figure 2 as functions of the applied potential drop (∆φ) measured with two Luggin electrodes, the positions of which are shown in Figure 1. The average DBL thickness estimated by Le´veˆque’s equation1,6 is about 63 µm. The limiting current density calculated by eq 4 is then 8.1 A m-2. The calculation of the total DBL thickness and its different zones near the MK-40 membrane starting from the experimental data presented in Figure 2 should be considered as an example of application of the developed method to describe ion-transfer phenomena at overlimiting current densities. 10.2. Calculation of the DBL Effective Thickness. To find the effective thickness of the DBL, δ′(u), it is sufficient to know + (u) the partial current densities of the Na+ (I(u) Na) and H (IH ) ions through the cation-exchange membrane. Remarking that the current density of the H+ ions through the membrane is equal to that of OH- ions through the DBL (I(u) OH), we can apply eq 36 under the form

I(u) Na )

2DNac(u) Na0F (u)

δ′

-

DNa (u) I DOH H

(44)

The results of calculation of δ′(u) as a function of I(u) according to eq 44, which is the modified Kharkats equation, are presented in Figure 4. As can be seen, δ′(u) decreases from δ(u) 0 ) 63 µm at the underlimiting currents to 12 µm attained at a current density nearly 10 times higher than the limiting one. This decrease may be attributed to coupled convection occurring as electroconvective mixing of the depleting solution near the interface.11-13,46,47 This mechanism is rather probable, taking into account the relatively high value of the SCR thickness (δSCR ) δ - δLEN) in comparison with the effective DBL thickness δ′ (Figure 4). The contribution of gravitational convection should be low in the considered system due to small concentration



Figure 4. Effective thickness of the DBL (δ′(u)) and the thickness of (u) ) as functions of the current density applied to the SCR (δSCR membrane system described in Figure 2.

gradients, small intermembrane spacing, small applied currents, and the relatively high velocity of solution flow.49 10.3. Calculation of the Thickness of the SCR and Its Different Zones. This calculation is made by fitting the total DBL thickness δ(u) C near the CEM when adjusting the theoretical, obtained by the integration of eqs 12 and 13, and experimental voltage values at a given current density. To calculate the potential drop (pd) between the tips of two measuring electrodes (∆φtot), one should evaluate the pd across the anion- (UA) and cation-exchange (UC) membranes, including the Donnan pd at the interfaces and the pd across the DBLs adjoining these membranes, as well as the pd (Usol) across the bulk solution layer between the external edges of the DBLs in the desalting compartment (Figure 1a)

∆φtot ) UA + UC + Usol

(45)

The two contributions from the Nernst electrode potentials cancel each other out. The resistance of the measuring circuit is very high, and the current flowing in this circuit is small. Hence, the concentration polarization of the electrodes is negligible, as well as the possible asymmetry of electrode potentials produced by the current flow. The calculation of the pd across the depleted DBL near the cation-exchange membrane was carried out after resolving the NPP eqs 12 and 13 under the QCD assumption. In this case, the dimensionless fluxes ji and the small parameter were found (u) from the experimental fluxes j(u) Na and jH at a given current (u) (u) density and a tentative value of δC , eq 3; j(u) Cl and jH were assumed to be zero in this DBL. The integration of the electric field found from eq 38 resulted in the pd between x ) 0 and x ) xs ) 1 - λ in the depleted DBL near the CEM. The pd over the depleted interface including the quasi-equilibrium layers in the solution and the membrane, from x ) 1 - λ to 1 + λh, was calculated from the Donnan eq 43. For the summary pd across both membrane interfaces (I and II), the following equation based on eq 43 was applied

∆φ(u) D

)

∆φID

+

∆φIID

II RT c1s ) ln I F c

(46)

1s

The boundary counterion concentration in the concentrating solution cII1s near the CEM was calculated by using the wellknown equation5,6,44

(

cII1s ) c10 1 +

I I IIlim

)

(47)

where IIIlim is the limiting current density determined by the DBL thickness in the concentrating compartment and c10 is the bulk counterion concentration in the concentrating compartment, being the same as that in the diluting one. The pd across the anion-exchange membrane, UA, with adjacent DBLs was not calculated by integrating eqs 12 and 13 as this calculation should require the thickness of the depleted DBL at this membrane, which involves, once more, fitting the parameter together with δ(u) C . Instead, UA was assumed to be proportional to UC, with a proportionality coefficient within the interval of [0.5, 1]. Effectively, UA should be lower than UC as the diffusion coefficient of Cl-, being the counterion for the AEM, is about 1.5 times higher than that of Na+, the counterion for the CEM. When assuming UA ) 0.5UC, eq 45 becomes

∆φ(u) tot

)

1.5U(u) C

(u) I (u)(h - δ(u) A - δC ) + κsol

(48)

where the second term in the right-hand part of eq 48 expresses U(u) sol, which is calculated by knowing the electric conductivity of the 0.002 M NaCl solution (κsol ≈ 0.025 S m-1) and applying Ohm’s law for the desalting compartment bulk. As the tips of the measuring electrodes were installed in the concentrating compartments at a small distance from the membranes, we have not included in ∆φ(u) tot the potential drops in these compartments, just the interfacial potential drops there, as given in eq 46. The ohmic potential drop inside of the CEM is negligible because of its relatively high conductivity in diluted solutions; the specific conductivity of the membrane in the 0.002 M NaCl solution is about 20 times higher than κsol. The value of δ(u) C is considered as true when the difference between the calculated and measured values of ∆φ(u) tot is (u) sufficiently small. After determining δC , the distribution of the electric field and concentrations can be established; then, the thickness of different zones of the DBL can be found. Figure 4 shows that at the highest current density applied in (u) the experiment (I(u) ) 9.52 mA cm-2, ∆φ(u) tot ) 7.7 V, UC ≈ 4 (u) (u) V, U(u) A ≈ 2 V, Usol ≈ 1.7 V), δSCR reaches 2.5 µm and makes up 17.5% of the overall thickness of the DBL (δ(u) C ) 14.2 µm) near the CEM. Note that if one assumes that the potential drop across the CEM (with the DBLs) is the same as that across the (u) (u) (u) AEM (U(u) C ≈ UA ≈ 3 V, Usol ≈ 1.7 V), δSCR will be equal to (u) 2.0 µm; δC ) 13.8 µm, whereas δ′(u) remains the same, 11.8 µm, at the same current density. Thus, the obtained estimation of the total DBL thickness (≈14 µm), where approximately 2.5 µm relates to the space charge region, seems sufficiently reliable for the given conditions of the membrane system, being a rather weak function of the assumption about the distribution of the potential drop. The contributions into the SCR are as follows: The transition layer makes up about 0.2 µm (found as 21/3), the migration zone constitutes 2.2 µm, and the quasi-equilibrium zone makes up 0.1 µm. These results show, in particular, a high importance of the current-induced convection (electroconvection as the most probable) in membrane systems, which is able to decrease the DBL thickness more than 5 times in the considered case. As for the thickness of the quasi-equilibrium zone (λ) that ranges between the point where c1 reaches its minimum and x ) 1, the numerical calculation for the conditions described above (I(u) ) 9.52 mA cm-2) gives λ(u) ≈ 0.11 µm. This value is close to that calculated by the Debye equation

14220 J. Phys. Chem. B, Vol. 111, No. 51, 2007 (u)

λ

)

x

RT(u) 0

2 2c(u) 1s (F)

Urtenov et al.

(49)

-6 M found from the numerical calculation with c(u) 1s ≈ 4.6 × 10 (u) λ ≈ 0.14 µm. Note that for the bulk concentration c(u) 10 ) 2 × 10-3 M, eq 49 yields λ(u) ≈ 0.007 µm.

11. Conclusion An approach allowing one to decouple the NPP equation system, in the case of steady-state 1D transport in a DBL adjacent to an electrode or an ion exchange membrane, is developed. The idea is to obtain equations for the electric field E not containing the concentration terms. The substitution of the E(x) function found from these equations into the NernstPlanck ones permits the concentration distribution to be found. This approach is rather general and suited for a multi-ionic system. In order to reduce the number of considered different ionic species, the grouping of the ions into valency classes is made. Then the supraconcentrations and suprafluxes are introduced into consideration as being the sum of products of concentrations (fluxes) and different powers of ionic charges. This results in obtaining recurrent relationships between (supra-) concentrations and fluxes, leading to an equation containing only one unknown variable, E(x). The cases of (pseudo) binary, ternary, and tertiary electrolyte systems are presented. Three differential equations, one for each electrolyte system, are deduced. This method, between others, gives an equation which links the ion fluxes at limiting and overlimiting currents. This equation, generalizing the known Kharkats equation, simplifies the analysis of the role of co-ions leaking through the membrane or generating in the course of a parallel chemical reaction. In particular, the influence of OH- ions produced by “water splitting” at a cation-exchange membrane and “exalting” the flux of salt cations toward the membrane surface is analyzed. The analysis of the NPP equations, in the case of the Na+, Cl-, and OH- ions present, has shown that at limiting and overlimiting currents, the DBL can be divided into four zones. In the first zone adjoining the bulk solution, the local electroneutrality is held. The distribution of the electric field and concentrations here can be obtained from the general equation system by setting the space charge density to be zero (the LEN assumption). The other zones form a space charge region near an IEM (or an electrode). In the second, migration zone, the charge density is no longer zero, but its variation there as well as the variation of concentrations is relatively small. The electroneutral and migration zones are separated by a transition zone. The transition zone arises at currents lower, but very close to, the limiting current when the migration zone is absent; the latter emerges only at overlimiting currents. The solution of the general equation system, the E(x) and ci(x) functions, can be found in the three indicated zones (electroneutral + transition + migration zones) by setting the derivative of the space charge density to be zero. This assumption, which implies a quasiuniform charge density distribution (the QCD assumption), may be considered as a new condition, allowing a more accurate, with respect to the LEN assumption, approximation for resolving the NPP equations. The QCD assumption is valid in the electroneutral and charged zones of the DBL, except for at the boundary quasi-equilibrium zone, which is a thin layer where the electric field and the concentration gradients are both high. E(x) and ci(x) can be approximately found here from the Nernst-Planck equations when the flux density term is neglected. Hence, the local equilibrium can be admitted, and the

potential drop there can be approximately calculated by the Donnan equation. The thickness of this zone is on the order of the Debye length determined by the minimum salt counterion concentration reached at the point separating the migration and quasi-equilibrium parts of the SCR. The thickness of all four zones is evaluated for a real membrane system as a function of the current density applied. Under the experimental conditions of a diluted solution (0.002 M NaCl), small intermembrane spacing (500 µm) and high average solution velocity (3.2 cm s-1), the thickness of the nondisturbed DBL is small (63 µm), and the parameter is relatively high (4.6 × 10-7). Under high concentration polarization (with the voltage across a cell pair of about 8 V), the effective thickness of the DBL, approximately equal to the electroneutral zone thickness, becomes about 12 µm. The thickness of the migration zone of the SCR is close to 2.2 µm; that of the transition zone, found as 21/3, is about 0.2 µm; the quasi-equilibrium zone remains thin, near 0.1 µm. Our model does not consider the origin of the DBL decrease under increasing voltage. It gives the means to evaluate the thickness of the DBL and its different zones by setting the model parameters in accordance with experimental current-voltage curves. The main attention in this study is paid to a membrane system under high polarization. However, it is of interest to consider the state where the current density is close to the limiting one. As the current-induced electroconvection partially destroys the DBL, the limiting current density increases with increasing current. Hence, the state under the limiting current is not a “point” state but rather a state related to an interval of current densities. The accurate consideration of this state requires careful analysis of the transition zone, in particular, the elucidation of the question about the applicability of the QCD assumption in this case. This and other questions remain for a following study. Acknowledgment. The research is supported by RFBR (Russian Foundation of Basic Research), Grants 06-03-96606-a and 07-08-00533-a, and RFTP (Russian Federal Targeted Programme), Contracts 02.513.11.3037 and 02.513.11.3163, as well as by INTAS, Grant 05-1000007-416. List of Symbols AbbreViations AEM, CEM, IEM: anion-, cation-, ion-exchange membrane, respectively DBL: diffusion boundary layer EDL: electric double layer LEN assumption: local electroneutrality assumption pd: potential drop QCD assumption: quasi-uniform distribution of the space charge density assumption SCR: space charge region Symbols a: integration constant, eq 17 c: electrolyte concentration ci: concentration of ion i Ci: group concentration, eq 10 di: dimensionless diffusion coefficient of ion i, eq 3 Di: diffusion coefficient of ion i E: electric field F: Faraday constant Gn,k: “supraflux”, eq 14 I: current density ji: flux density of ion i

Decoupling of the Nernst-Planck and Poisson Equations Ji: group flux density, eq 11 n: number of groups, eq 10 ni: number of ions in group i, eq 10 N: total number of ions LD: Debye length, eq 3 qn,k: parameter, eq 15 R: universal gas constant Sn,k: “supraconcentration”, eq 14 ti: transport number of ion i in solution T: absolute temperature Ti: effective transport number of ion i in the membrane, eq 30 UA, UC, Usol: pd across anion- and cation-exchange membrane, and solution bulk, respectively x: coordinate normal to the membrane zi: charge number of ion i Greek Symbols δ: total diffusion layer thickness δ′: effective thickness of the diffusion layer δLEN: thickness of the electroneutral zone of the diffusion layer δSCR: thickness of the space charge region in the diffusion layer (u): dielectric permittivity, eq 2 : small parameter, eq 3 κ: electric conductivity λ: thickness of the quasi-equilibrium space charge zone in solution λh: thickness of the quasi-equilibrium space charge zone in the membrane λtr: half of the thickness of the transition space charge zone in solution F: electric space charge density η,ξ: variables, eq 26 φ: electric potential ∆φ: potential difference Indices Subscripts i: ion number lim: related to the limiting current 0: related to the solution bulk s: related to the “surface” equal to the point where the salt counterion concentration attains its minimum 1: counterion 2: co-ion Superscripts I and II: left-hand and right-hand interface, respectively (u): dimension quantity (having a unit) References and Notes (1) Newman, J. S. Electrochemical Systems; Prentice Hall: Englewood Cliffs, NJ, 1973. (2) Lakshminarayanaiah, N. Equations of Membrane Biophysics; Academic: New York, 1984. (3) McKelvey, J. P. Solid State and Semiconductor Physics; Krieger: Malabar, FL, 1982. (4) Timashev, S. F. Physical Chemistry of Membrane Processes; Ellis Horwood Ltd.: London, 1990. (5) Rubinstein, I. Electro-diffusion of Ions; SIAM: Philadelphia, PA, 1990; p 254. (6) Zabolotsky, V. I.; Nikonenko, V. V. Ion Transport in Membranes (in Russian); Nauka: Moscow, 1996; p 390. (7) Lebedev, K.; Mafe´, S.; Alcaraz, A.; Ramirez, P. Chem. Phys. Lett. 2000, 326, 87. (8) Mafe´, S.; Ramirez, P.; Alcaraz, A.; Aguilella, V. Handbook on Bipolar Membrane Technology; Kemperman, A. J. B, Ed.; Twente University Press: Enschede, The Netherlands, 2000; p 49.

J. Phys. Chem. B, Vol. 111, No. 51, 2007 14221 (9) Rubinstein, I.; Shtilman, L. J. Chem. Soc., Faraday Trans. 1979, 75, 231. (10) Dukhin, S. S. AdV. Colloid Interface Sci. 1991, 35, 173. (11) Mishchuk, N. A.; Takhistov, P. V. Colloids Surf., A 1995, 95, 119. (12) Rubinstein, I.; Zaltzman, B. Phys. ReV. E 2000, 62, 2238. (13) Rubinstein, I.; Zaltzman, B.; Lerman, I. Phys. ReV. E 2005, 72, 011505/1. (14) Wang, S.-C.; Lai, Y.-W.; Ben, Y.; Chang, H.-C. Ind. Eng. Chem. Res. 2004, 43, 2902. (15) Urbanski, J. P.; Thorsen, T.; Levitan, J. A.; Bazant, M. Z. Appl. Phys. Lett. 2006, 89, 143508. (16) Barany, S.; Mishchuk, N. A.; Prieve, D. C. J. Colloid Interface Sci. 1998, 207, 240. (17) Ben, Y.; Demekhin, E. A.; Chang, H.-C. J. Colloid Interface Sci. 2004, 276, 483. (18) Guzma´n-Garcia, A. G.; Pintauro, P. N.; Verbrugge, M. W.; Hill, R. F. AIChE J. 1990, 36, 1061. (19) Cervera, J.; Manzanares, J. A.; Mafe´, S. J. Membr. Sci. 2001, 191, 179. (20) Ariel, N.; Ceder, G.; Sadoway, D. R.; Fitzgerald, E. A. J. Appl. Phys. 2005, 98, 023516. (21) Bazant, M. Z.; Chu, K. T.; Bayly, B. J. SIAM J. Appl. Math. 2005, 65, 1463. (22) Levich, V. G. Dokl. Akad. Nauk SSSR 1949, 67, 307. (23) Levich, V.G. Physicochemical Hydrodynamics; Prentice Hall: Englewood Cliffs, NJ, 1962. (24) Grafov, B. M.; Chernenko, A. A. Zh. Fiz. Khim. 1963, 37, 664. (25) Newman, J. J. Chem. Soc., Faraday Trans. 1965, 61, 2229. (26) MacGillivray, A. D. J. Chem. Phys. 1970, 52, 3126. (27) Dukhin, S.S.; Deryagin, B.V. Electrophoresis (in Russian); Mir: Moscow, 1976. (28) Smyrl, W. H.; Newman, J. Trans. Faraday Soc. 1967, 63, 207. (29) Rubinstein, I.; Zaltzman, B. Math. Models Methods Appl. Sci. 2001, 11, 263. (30) Zaltzman, B.; Rubinstein, I. J. Fluid Mech. 2007, 579, 173. (31) Mafe´, S.; Pellicer, J.; Aguillela, V. M. J. Comput. Phys. 1988, 75, 1. (32) Manzanares, J. A.; Murphy, W. D.; Mafe´, S.; Reiss, H. J. Phys. Chem. 1993, 97, 8524. (33) Chu, K. T.; Bazant, M. Z. SIAM J. Appl. Math. 2005, 65, 1485. (34) Zabolotskii, V. I.; Nikonenko, V. V.; Korzhenko, N. M.; Seidov, R. R.; Urtenov, M. Kh. Russ. J. Electrochem. 2002, 38, 810. (35) Zabolotskii, V. I.; Lebedev, K. A.; Lovtsov, E. G. Russ. J. Electrochem. 2006, 42, 836. (36) Listovnichy, A. V. Elektrokhimiya 1989, 25, 1651. (37) Nikonenko, V. V.; Zabolotskii, V. I.; Gnusin, N. P. Elektrokhimiya 1989, 25, 301. (38) Urtenov, M. Kh.; Nikonenko, V. V. Elektrokhimiya 1993, 29, 239. (39) Nikonenko, V. V.; Urtenov, M. Kh. Elektrokhimiya 1996, 32, 195. (40) Babeshko, V. A.; Zabolotskii, V. I.; Korzhenko, N. M.; Seidov, R. R.; Urtenov, M. Kh. Russ. J. Electrochem. 1997, 33, 793. (41) Chu, K. T. Asymptotic Analysis of Extreme Electrochemical Transport. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, 2005. (42) Babeshko, V. A.; Zabolotskii, V. I.; Kirillova, E. V.; Urtenov, Kh, M. Dokl. Ross. Akad. Nauk 1995, 344, 485. (43) Babeshko, V. A.; Zabolotskii, V. I.; Seidov, R. R.; Urtenov, M. Kh. Russ. J. Electrochem. 1997, 33, 785. (44) Manzanares, J.; Kontturi, K. In Encyclopedia of Electrochemistry, Vol. 2, Interfacial Kinetics and Mass Transport, Diffusion and Migration; Bard, A. J., Stratmann, M., Calvo, E. J., Eds.; Whiley Publishing Inc.: Indianapolis, IN, 2003; pp 81-121. (45) Zabolotsky, V. I.; Nikonenko, V. V.; Pismenskaya, N. D.; Laktionov, E. V.; Urtenov, M. Kh.; Strathmann, H.; Wessling, M.; Koops, G. H. Sep. Purif. Technol. 1998, 14, 255. (46) Belova, E. I.; Lopatkova, G. Yu.; Pismenskaya, N. D.; Nikonenko, V. V.; Larchet, C.; Pourcelly, G. J. Phys. Chem. B 2006, 110, 13458. (47) Pismenskaya, N. D.; Nikonenko, V. V.; Belova, E. I.; Lopatkova, G. Yu.; Sistat, Ph.; Pourcelly, G.; Larshet, K. Russ. J. Electrochem. 2007, 43, 307. (48) Balster, J.; Yildirim, M. H.; Stamatialis, D. F.; Ibanez, R.; Lammertink, R. G. H.; Jordan, V.; Wessling, M. J. Phys. Chem. B 2007, 111, 2152. (49) Zabolotsky, V. I.; Nikonenko, V. V.; Pismenskaya, N. D. J. Membr. Sci. 1996, 119, 171. (50) Shaposhnik, V. A.; Vasil’eva, V. I.; Ugryumov, R. B.; Kozhevnikov, M. S. Russ. J. Electrochem. 2006, 42, 531. (51) Simons, R. Electrochim. Acta 1985, 30, 275. (52) Ibanez, R.; Stamatialis, D. F.; Wessling, M. J. Membr. Sci. 2004, 239, 119. (53) Umnov, V. V.; Sheldeshov, N. V.; Zabolotskii, V. I. Russ. J. Electrochem. 1999, 35, 411. (54) Kilic, M. S.; Bazant, M. Z. Phys. ReV. E 2007, 75, 021503.

14222 J. Phys. Chem. B, Vol. 111, No. 51, 2007 (55) Vasil’eva, A. B.; Butuzov, V. F. Asymptotic Decomposition of Solutions for Singularly Perturbed Equations (in Russian); Nauka: Moscow, 1973. (56) Vasil’eva, A. B.; Butuzov, V. F.; Kalachev, L. V. The Boundary Function Method for Singular Perturbation Problems; SIAM: Philadelphia, PA, 1995. (57) Boglaev, J. P. Zh. Vyssh. Mat. 1970, 10, 958. (58) Schlo¨gl, R. Z. Phys. Chem. (Frankfurt/Main, Ger.) 1954, 1, 305. (59) Kramers, H. A. Physica 1940, 7, 284. (60) Oldham, K. B.; Feldberg, S. W. J. Phys. Chem. B 1999, 103, 1699. (61) Urtenov, M. Kh. Boundary-Value Problems for the NernstPlanck-Poisson Equations (in Russian); Kuban State University: Krasnodar, Russia, 1988 and 1999. (62) (a) Collatz, L. Functional Analysis and Numerical Mathematics; Academic Press Inc.: New York, 1966. (b) Ortega, J. M.; Rheinboldt, W.

Urtenov et al. C. IteratiVe Solution of Nonlinear Equations in SeVeral Variables; Academic Press: New York, 1987. (63) Buck, R. P. J. Membr. Sci. 1984, 17, 1. (64) Heyrovsky, J.; Kuta, J. Principles of Polarography; Academic Press: New York, 1966. (65) Damaskin, B. B.; Petriy, O. A. Introduction to the Electrochemical Kinetics (in Rusian); Vysch. Shkola: Moscow, 1975, p 173. (66) Kharkats, Yu. I. SoV. Electrochem. 1985, 21, 917. (67) Kharkats, Yu. I.; Sokirko, A. V. J. Electroanal. Chem. 1991, 303, 27. (68) Rubinstein, I. React. Polym. 1984, 2, 117. (69) Strathmann, H. Ion-Exchange Membrane Separation Processes; Membrane Science and Technology Series 9; Elsevier: Amsterdam, The Netherlands, 2004, p 348.

Decoupling of the NernstâPlanck and Poisson Equations. Application

Recommend Documents

Decoupling of the NernstâPlanck and Poisson Equations. Application