## Abstract

The electric potential of counter-ions (protons) in an infinite cylindrical channel is presented as a solution of the Poisson–Boltzmann equation, involving a constant ion charge density along the wall. The distribution of protons is derived and used subsequently to compute the velocity profile and mass flow rate of the corresponding electro-osmotic flow, driven by an electric field. Analytical expressions are derived for all quantities, including the conductivity and water drag coefficient. This analysis relates especially to cylindrical nano-channels of polymer electrolyte membranes such as Nafion and addresses the validity of continuum models for these materials.

## 1. Introduction

Following the work of Gebel (2000), Schmidt-Rohr & Chen (2008) recently proposed that the micelle structure of polymer electrolyte membranes (PEMs) might consist of only cylindrical nano-channels, which facilitate water and proton transport, rather than large water pore clusters connected by smaller nano-channels as in Gierke’s model (Thompson *et al.* 2006). This raises the question as to whether, and how, water and proton transport in PEMs can be modelled by electro-osmotic flow in cylindrical pores (Eikerling & Kornyshev 2001). In these materials, negatively-charged ions are located at the channel walls. These form an ion-exchange equilibrium with the surrounding water molecules, resulting in a mixture of bound and dissolved protons.

Owing to the small diameter of the channels (about 3 nm), it is unclear whether continuum models or Lattice–Boltzmann simulations are required for such flows (Akinaga *et al.* 2008). Qiao & Aluru (2003) compared both approaches and concluded that continuum models, e.g. the Poisson–Boltzmann (PB) equation coupled to Stokes flow, might indeed be able to describe such systems, provided that the flow near the channel wall is modelled correctly. This is consistent with the conclusions of Eikerling & Kornyshev (2001) who predict a layer near the wall, consisting of rather tightly bound protons, and an interior layer of more mobile protons. However, strictly speaking, the PB equation describes only stationary electro-diffusion systems, meaning the absence of a non-convective ion flux (Rubinstein 1990). Hence, the validity of such equilibrium approaches to non-equilibrium settings needs to be examined carefully.

The PB equation has received much attention in the literature, as it describes many electrochemical systems or electrostatics near charged surfaces (Rubinstein 1986, 1990; Van Aken *et al.* 1990; Borukhov *et al.* 2000; Ohshima 2003; Dobrynin & Rubinstein 2005). This applies especially to rectangular and cylindrical geometries (Fouss *et al.* 1951; Hinderberger *et al.* 2005; Chen *et al.* 2008). Usually, it requires numerical techniques to solve the equations (Qiao & Aluru 2004), whereas analytical, meaning closed form, solutions can rarely be found (Rice & Whitehead 1965; Tsao 1998; Lipkowitz *et al.* 2003). For this reason, the PB equation is often linearized (Tseng *et al.* 2005). For a review regarding the PB equation, the reader may be referred to Grochowski & Trylska (2007) and Lipkowitz *et al.* (2003).

Following the ideas of Qiao & Aluru (2003) related to a rectangular two-dimensional channel (i.e. with finite height but with infinite width), this paper expands a model for electro-osmotic flow to a three-dimensional cylindrical domain, as shown in figure 1. The channel walls have a uniform surface charge *σ*_{s}, balanced by a proton concentration *c* inside the channel so that global electro-neutrality is preserved. This serves as a simplified model for a nano-pore in polymer electrolyte membranes such as Nafion, where the surface charges are partially ionized sulphonic acid groups (SO_{3}^{−}), subject to an ion-exchange equilibrium. Eikerling *et al.* (2001) refer to this concept as a ‘proton cloud’, in which the rapid formation of various proton–water complexes and the transition between them are taking place so fast that an effective PB equation is considered sufficient to describe the dynamics. However, although in their paper an effective proton mobility was computed based on transition probabilities, we will pursue the concept of Stokes flow in a Nafion nano-pore.

Owing to the symmetry and the constant pore radius of the geometry in this paper, the PB equation might indeed be a good *quasi-equilibrium* approximation for the full electro-osmotic flow problem (Rubinstein 1990). However, when discrete wall charges are considered (Eikerling & Kornyshev 2001) or the channel contains a variation in radius along its lateral *z*-axis, the PB equation could be a poor approximation. This depends on how the timescale for convective flow down the channel compares with that of the diffusive radial flow of ions which seeks to establish locally the thermodynamic equilibrium assumed in the derivation of the PB equation. With these potential limitations in mind, we will proceed in our analysis.

Regarding the PB equation, the model presented in this paper has been solved, in part, previously by Tsao (1998), who used a perturbation ansatz, Lipkowitz *et al.* (2003) and Eikerling *et al.* (2001). However, we propose a new method for the derivation of the solution, leaning on the ideas of Fouss *et al.* (1951) who solved analytically for the potential between two concentric cylinders. Moreover, to the best of the authors’ knowledge, the analytical solutions of the corresponding Stokes flow problem, including the velocity profile, mass flow rate, conductivity and water drag, are derived for the first time. This is the main point of this paper.

Sections 2 and 3 will present the system of equations, namely PB equation and Stokes flow, and their analytical solution, containing the velocity profile, mass flow rate, conductivity and water drag. We conclude our analysis in §4.

## 2. The model

We study an infinite cylindrical pore with a uniform surface charge density *σ*_{s}. We assume that the flow is driven in the *z*-direction by only an external field, not by a pressure gradient. Hence, there is no dependency of the solutions for the electric potential *ψ*, the proton concentration *c* or the velocity along the channel *u* on the angular or lateral coordinates (*α*,*z*).

It is important to mention that unless a uniform surface charge density is assumed, as well as global electro-neutrality, the electric field does not vanish outside the pore, and the problem needs to be solved in infinite space (Eikerling & Kornyshev 2001).

### (a) The governing equations

Following Qiao & Aluru (2003), we have to solve the following system of equations:
2.1
and
2.2
with the corresponding boundary conditions (BCs)
2.3
2.4
2.5
and
2.6
Here, the vacuum permittivity has the value *ε*_{0}=8.854×10^{−12} C^{2}/(N m^{2}), *q* is the electron charge (1.6×10^{−19} C), is the valency of the protons (), *c*_{0}=*c*(0) is the proton concentration at the channel centre, *k*_{B} is the Boltzmann constant, *T* is the temperature, *μ* is the viscosity, *E*_{ext} is the (external) electric field and *σ*_{s} is the charge density along the wall. In addition, *r*=0 denotes the channel centre and *r*=*R* denotes the channel wall. Note that we neglected a variation of the viscosity (*μ*=3.35×10^{−4} Pa s) or dielectric constant, *ε*, with *r*.

By applying the PB equation (2.1), we assume that (i) the diffusive flux balances the migration flux in the Nernst–Planck equation if one was to apply the complete set of Poisson–Nernst–Planck (PNP) and Navier–Stokes (NS) equations (Rubinstein 1990) and consequently (ii) convection dominates the proton flow. Moreover, we assume at first that the pressure only plays a minor role in the movement of water along the channel and, hence, can be neglected. Therefore, it does not appear in the Stokes flow equation (2.2). We will come back to this assumption later in this paper. It should also be said that the incompressibility condition for water and the conservation of protons is fulfilled because of the translational invariance of the domain, being an infinite straight cylinder.

The reader should be reminded that in using the PB equation, the proton concentration is determined thermodynamically by 2.7

The BCs exploit the symmetry of the system. Equations (2.3) and (2.4) are derived from Gauss’ law, using the divergence theorem, whereas equations (2.5) and (2.6) are simply a no-slip and a symmetry condition, respectively. The relation (2.4) holds only if global electro-neutrality is assumed, resulting in a vanishing electric field for *r* > *R*. The balance between total proton and ion charges in a cylindrical segment of length *L* reads
2.8
This can be shown to be equivalent to the BC (2.4) by, first, integrating the left-hand side of equation (2.8) via substitution of the left-hand side of equation (2.1) and, second, by utilizing the divergence theorem and the symmetry of the problem.

We will solve the problem in two steps. First, the PB equation (2.1) is solved, which also yields the proton concentration via the above equation (2.7) and, subsequently, the Stokes flow (2.2) is computed.

### (b) A comment on gauge invariance

Using the Laplace operator in cylindrical coordinates, the PB equation under consideration is (Tsao 1998)
2.9
The boundary conditions take the form (prime denotes the first derivative and double prime the second derivative)
2.10
and
2.11
In addition, we choose
2.12
This third condition does not lead to an over-determined problem. In fact, the problem posed is gauge invariant and this choice is possible because the right-hand side of equation (2.9) can always be re-written according to
2.13
with
2.14
Hence, we can always re-scale *c*_{0} and choose *ψ*(0)=0.

Another way of interpreting the situation is to consider *c*_{0} to be an eigenvalue of equation (2.9), subject to the conditions (2.10)–(2.12).

## 3. Solution

### (a) PB equation

We will now solve equation (2.9), subject to the BCs (2.10)–(2.12). This problem was solved previously by Tsao (1998) with a power series approach and by Lipkowitz *et al.* (2003), who make use of the corresponding solution for two concentric cylinders. However, we will pursue a different approach by using direct integration, following the ideas of Fouss *et al.* (1951).

In order to non-dimensionalize the equations, we first set
3.1
Then, we find
3.2
and
3.3
The second step in the non-dimensionalization procedure is the variable transformation *r*=*x**R*, where *R* is the radius of the pore. We obtain
3.4
where we defined a dimensionless parameter
3.5
which relates to the Debye length
3.6
The BCs now read
3.7
We re-write the differential equation
3.8
and set (Fouss *et al.* 1951)
3.9
(Note that special attention is required at *x*=0.) This yields
3.10
and upon substitution into the differential equation
3.11
with the new BCs
3.12
3.13
and
3.14
Next, we set and *ϱ*(*y*(*x*))=*ϕ*(*x*). Then *ϱ* satisfies
3.15
and, consequently, the three BCs for *ϕ* transform into
3.16
where we used the scaling of *ϕ*(*x*) and as *x*→0 (see equations (3.12) and (3.13)). Substituting equation (3.15) into equation (3.11) yields
3.17
Finally, we choose *θ*=*ϱ*+*y* to obtain
3.18
with the BCs
3.19
3.20
and
3.21
Integration and use of the BCs give
3.22
where
3.23
As we know that *θ*′>0, we have
3.24
Let *θ*_{0} denote the value of *θ* when *y*=0. Our BC (3.21) now reads
3.25
Integration of equation (3.24) gives
3.26
Solving for *θ*(*y*) and using equation (3.25), we find
3.27
Making all substitutions in reverse order, we obtain
3.28
The last step is now to ensure that all three BCs in equation (3.7) are fulfilled. Imposing and Taylor expanding about *x*=0, we find
3.29
Hence, the argument must be 1. In other words, we need
3.30
from which we see that
3.31
Substituting the expression for *λ*, we obtain
3.32
We discover the obvious results that *σ*_{s}(0)=0 and . However, we also recover the well-known phenomena of maximum charge density according to
3.33
Using the relation (3.32) in equation (3.28), our solution assumes the simple form
3.34
We see immediately that *λ* is restricted to the range 0≤*λ*<8. This solution was also found by Tsao (1998), Eikerling *et al.* (2001) and (apart from a sign error) Lipkowitz *et al.* (2003) but not with direct integration.

As and hold, we also know that another BC, namely , is fulfilled because is differentiable at *x*=0. This can be verified by computing the derivative explicitly from the solution (3.34).

This only leaves the last boundary condition 3.35 It can again be verified easily that it is also fulfilled, using the relation (3.31) 3.36

It should be mentioned here that although the solution (3.34) is finite at *x*=0, corresponding to *r*=0, the work of Fouss *et al.* (1951) for a cylindrical charged molecule contains a singularity at *r*=0 which lies outside the physical domain. However, the techniques involved in each approach make use of the same transformations.

Lastly, the proton concentration can be computed by substituting equation (3.34) into the non-dimensionalized analogue of equation (2.7), yielding
3.37
A plot of *c*(*x*) is provided in figure 2 for different values of *λ*.

### (b) The reference case

As in most electro-osmotic flow problems, a natural lengthscale, the Debye length *d* in equation (3.6), emerges. For a typical water-filled pore in a PEM such as Nafion, we can estimate this lengthscale. We first solve equation (3.32) for *c*_{0}:
3.38
and use the following estimates for the remaining parameters:

—

*R*=10^{−9}m (Schmidt-Rohr & Chen 2008);—

*T*=353 K (operating temperature of PEM fuel cells);—

*σ*_{s}=−*q*/(1.0×10^{−9}m)^{2}(Eikerling & Kornyshev 2001);—

*ϵ*=45 (Paddison*et al.*1998), based on the rather high water/sulphonic-acid-group ratio of about 10 which is given by*R*and*σ*_{s}; and— V m

^{−1}, based on a voltage drop of 0.1 V across a Nafion membrane of thickness 50*μ*m.

Using the above expressions in equation (3.6), we find a Debye length of
3.39
which is larger than the value found by Eikerling & Kornyshev (2001). However, this lengthscale might change slightly with the use of a non-constant permittivity (Qiao & Aluru 2003)
3.40
or discrete surface charges (Eikerling & Kornyshev 2001). On the other hand, it seems to reconfirm the validity of the PB equation following the work of Corry *et al.* (2000). The authors found that PNP or PB model predictions resemble those of Brownian dynamics when the channel radius exceeds two Debye lengths.

Looking again at figure 2, the value of *λ*=6.13 corresponds to a situation where the protons are overwhelmingly concentrated near the channel wall and, therefore, near the sulphonic acid groups. This will have important consequences for the model results below.

### (c) Stokes flow

Finally, we are able to integrate the Stokes flow equation (2.2) analytically, using the solution (3.37). In non-dimensional cylindrical coordinates (*r*→*x**R*), we need to solve
3.41
where we assume that there is only a *z*-component *u*(*x*) of the velocity that solely depends on the radial coordinate *x*. We substitute equation (3.37) for *c*(*x*) and re-write the equation as
3.42
Integration on both sides and subsequent division by *x* results in
3.43
The constant of integration is determined by the symmetry condition (d*u*/d*x*)(0)=0. In order for the right-hand side of equation (3.43) to remain finite at *x*=0, we need *C*_{1}=4*c*_{0}*q**E*_{ext}*R*^{2}/*λ**μ*. This leads to
3.44
where the right-hand side does indeed tend to zero as *x*→0. We integrate again and find
3.45
With the no-slip condition at the channel wall, *u*(1)=0, the integration constant is found to be
3.46
Upon substitution, the velocity profile reads
3.47
from which we find the velocity at the channel centre
3.48
Three conclusions can be drawn. First, the velocity profile resembles that of Poiseuille flow when *λ* is small and we can Taylor expand the logarithm in equation (3.47) as
3.49
This approximation is valid for values of *λ* less than 0.5, as shown in figure 3.

Second, the average water velocity for our reference case is
3.50
This would mean that it takes about 0.4 s for an average water molecule to be dragged across a Nafion membrane of thickness 50 *μ*m if it consisted of only such cylindrical pores. This is smaller than diffusion timescales of about 20 s for water in Nafion with a thickness of 50 *μ*m that have been measured experimentally (Berg *et al.* 2004). However, we are dealing here with two different transport mechanisms, namely water drag versus water diffusion. In the context of our PB–Stokes flow model, the latter can only be modelled by a pressure gradient (Berg *et al.* 2004) and this will be addressed further below.

Third, the total mass flux, , through a cross section of a cylinder can be computed according to
3.51
with *ρ* being the density of the water. The substitution of equation (3.47) into the integral yields an analytical expression for :
3.52
where *f*(*λ*) represents the expression in the brackets of the second line in equation (3.52). This has a (Poiseuille) limit for small *λ* of
3.53
We observe for small, permissible values of *λ* that the mass flux always scales like *R*^{4}, which is the same as for the classical Poiseuille flow. However, the *R*^{4} scaling is not valid for larger values of *λ*.

### (d) Conductivity and water drag

To gauge the validity of the continuum approach, the model predictions need to be compared to experimental results such as those for conductivity and water drag. This will be investigated next.

The total proton flux, or current, is determined by
3.54
having unit ampere. We can write the right-hand side as
3.55
similar to the expression for , where *g*(*λ*) consists of an integral that depends on *λ* only
3.56
where we used integration by parts to solve the integral. In the limit *λ*→0, we obtain *g*(0)=1/32.

We can now compute two quantities. First of all, the conductivity *σ*_{p} is determined by
3.57
and for our reference case it is found to be *σ*_{p}=0.15 S cm^{−1}. It appears too large by a factor of 2–3 when compared with literature values (Paddison *et al.* 1998). However, one should really divide the right-hand side of equation (3.57) by the effective cross section that each pore occupies which is easily two times as large (Schmidt-Rohr & Chen 2008) when the hydrophobic backbone is taken into account. This brings the effective conductivity in line with experimental data (*σ*_{p}=0.06 S cm^{−1}).

The second quantity of interest is the electro-osmotic drag, *n*_{drag}. In this model, it is given by the ratio of water flux to proton flux. Dividing the respective expressions (3.52) and (3.55), taking into account that we need to convert both expressions to molar fluxes, yields
3.58
Here, we introduced the molar weight of water, *w*=0.018 kg mol^{−1}. Figure 4 shows the water drag as a function of *λ*, given by equation (3.58). Note that our reference case corresponds to *n*_{drag}=22.7, i.e. an order of magnitude higher than literature values (Zawodzinski *et al.* 1995). This obviously stems from the fact that the protons are concentrated near the walls, where the fluid velocity is low, whereas their concentration at the channel centre is small, where *u* is large. Note that the proton concentration at the centre of the pore, *c*_{0}, generally decreases with an increasing wall charge *σ*_{s}, resulting in a smaller *λ*. Hence, we find larger water drags at smaller values for *λ*. In the limit *λ*→0, the water drag is determined by *n*_{drag}=*ρ**F*/*q**c*_{0}*w*=71.9.

Another point is that we neglected pressure and also assumed that the proton flux will be dominated by convection and not migration. Therefore, one could try and remedy this water-drag discrepancy through a reduction of the fluid velocity at the centre of the pore by introducing a pressure gradient along the pore. This is shown in figure 3 where a constant pressure gradient *δ**P* was introduced in the Stokes equation (3.41), still allowing for an analytical expression for *u*(*x*). However, this requires unrealistically large pressure gradients so as to have an impact on *n*_{drag}, equivalent to approximately 100 atm. across the Nafion membrane. On the other hand, applying a pressure gradient equivalent to 1 atm. across a 50 *μ*m Nafion membrane, neglecting any external electric fields, results in an average channel velocity of about 5×10^{−7} m s^{−1}. This means that an average pressure-driven water molecule would cross the membrane in about 100 s, which is the correct order of magnitude for water diffusion across a PEM (Berg *et al.* 2004).

The other, more successful approach is to compute the total proton flux as the sum of convective flux and migration, following the Nernst–Planck equation, according to
3.59
(Note that diffusion does not play a role as the concentration gradient of protons along the pore is zero to leading order.) In this case, we find a water drag of *n*_{drag}=1.26 when we choose the proton diffusion coefficient to be *D*=10^{−8} m s^{−1}. This value seems to be too high by an order of magnitude though (Paddison *et al.* 2000). Therefore, the issue of excessively large water drags in our simple continuum model has got no easy fix and needs to be investigated in greater detail in future research. If this discrepancy cannot be resolved, there is a chance that continuum models are incapable of describing the water and proton transport in PEM nano-pores.

One explanation for this discrepancy could be that continuum models might overestimate the water velocity at the centre (see fig. 6 in Qiao & Aluru 2003). Another one could be that the distribution of sulphonic acid groups needs to be treated in a discrete manner (Eikerling & Kornyshev 2001).

### (e) Sensitivity analysis

The last question to be addressed is how a variation of parameters in the model, such as pore radius, dielectric constant or surface charge, affect the results. One parameter that we do not have to check is the external electric field because it cancels out in the derivation of the conductivity and the water drag.

We find that when the radius is varied between 0.5 and 1.5 nm (Schmidt-Rohr & Chen 2008), the surface charge between −*q*/(0.7×10^{−9} nm)^{2} and −*q*/(1.2×10^{−9} nm)^{2} (Eikerling & Kornyshev 2001) and the dielectric constant between 10 and 55 (Paddison *et al.* 1998), *λ* varies between 3.9 and 7.8. The lower value of 3.9 corresponds to the limit of continuum modelling (Corry *et al.* 2000).

This results in a conductivity that increases or decreases by a factor of 2 or 0.5, respectively. In other words, it remains close to literature values when the effective cross-sectional area of a pore is taken into account, as mentioned above.

Regarding the water drag, which is plotted in figure 4, the smallest value is achieved for *λ*=7.8 and it is about 5 and, hence, still too large. However, this includes only proton transport by convection and needs to be further investigated in a model that contains the full PNP equations.

## 4. Discussion and conclusion

This paper examines an electro-osmotic flow problem in an infinite cylindrical pore with a uniform surface charge density. It employs the PB equation and Stokes flow. Assuming constant permittivity, viscosity and radius of the pore, the governing equations can be solved in closed form. This reveals the distribution of the electric potential and the counter-ions (protons), the velocity profile of the water flow and its associated total flux, as well as the protonic current, conductivity and water drag.

The novelty of this contribution lies in the derivation of fluxes, conductivity and water drag. To the best of the authors’ knowledge, direct integration is employed for the first time, which differs in its approach from previous work. It results in analytical expressions for all quantities.

It is found that the electro-osmotic flow resembles that of the classical Poiseuille flow if (and only if) the dimensionless parameter *λ*, as defined in equation (3.5), is small. This is equivalent to a narrow pore.

The question of the applicability of the analysis in regards to polymer electrolyte membranes arises. First, the PB equation is, strictly speaking, valid only in systems with vanishing diffusion–migration, where convection dominates the ion flux (Rubinstein 1990). However, the PB equation is generally a popular approach to investigate non-stationary electro-osmotic systems (Eikerling & Kornyshev 2001; Qiao & Aluru 2003). For an infinite cylindrical pore with constant radius and constant, uniform surface charge density, the PB equation might indeed be applicable due to the translational invariance of the governing equations (*z*-direction).

However, a realistic PEM nano-channel is finite, varies in diameter along its main axis and has a discrete distribution of charges (Eikerling & Kornyshev 2001). This breaks the translational invariance of the governing equations. Therefore, it is crucial to assess how far the system deviates from equilibrium when a flux is present. One needs to compare the timescales of convective flow and diffusive proton flux. If the protons are given sufficient time to adjust locally to the thermodynamic equilibrium prescribed by equation (2.7), equivalent to the diffusive timescale being smaller than the convective timescale, the PB equation might provide a good quasi-equilibrium approximation for the full, non-stationary governing equations (Rubinstein 1990). According to Corry *et al.* (2000), a prerequisite for the validity of continuum equations in such channels is that the Debye length is less than half the channel radius. This is the case in our study that describes polymer electrolyte membranes such as Nafion.

At this point, our analysis focuses on the predictions of a simplified model, based on the PB equation coupled to Stokes flow. It is found that the only significant discrepancy between theoretical and numerical results lies in the value of *n*_{drag}. It indicates that our simple approach might be insufficient and that a more comprehensive model is needed, based on the full PNP equations (coupled to Stokes flow).

Therefore, future work will investigate the issue of large water drags as they occur in this paper. A radial-dependent viscosity and/or dielectric constant might be required for the model (Gupta *et al.* 1995; Hietala *et al.* 1997), along with a discrete distribution of sulphonic acid groups. In addition, the proton–wall and water–wall interaction might need to be modelled with greater care (Qiao & Aluru 2003). This would prevent the derivation of analytical expressions and relations, and call for numerical solution techniques. Last but not least, there exists the possibility that continuum models inherently fail to describe the proton dynamics as the latter are controlled by molecular rearrangements in a hydrogen-bonded network.

## Acknowledgements

The authors thank I. Rubinstein for the discussions and Toyota Motor Corporation for the financial support of this research.

## Footnotes

- Received February 6, 2009.
- Accepted May 19, 2009.

- © 2009 The Royal Society