Downloaded from http://rspa.royalsocietypublishing.

org/ on June 22, 2015

Proc. R. Soc. A (2011) 467, 3157–3169

doi:10.1098/rspa.2011.0080
Published online 15 June 2011

Analytical solution of the

Poisson–Nernst–Planck–Stokes equations
in a cylindrical channel
BY PETER BERG* AND JUSTIN FINDLAY
Faculty of Science, University of Ontario Institute of Technology, Oshawa,
Ontario, Canada L1H 7K4

We present a complete analytical solution for electro-kinetic flow of an acidic solution

in an infinite, circular, cylindrical channel. The flow of protons and water is described
by Poisson–Nernst–Planck equations, coupled to Stokes flow, while negative charges,
situated along the walls, maintain electro-neutrality. Several system characteristics are
computed in closed form such as the velocity profile, mass flux, current and water drag.
The two latter quantities allow for a comparison of analytical results with experimental
data of systems such as polymer electrolyte membranes.
Keywords: Poisson–Nernst–Planck equations; Stokes flow; analytical solution;
cylindrical channel; polymer electrolyte membranes

1. Introduction

Nanofluidics plays an increasingly important role in the research and development

of new technologies. Examples are novel electrode and electrolyte materials for
fuel cells and batteries. In this context, the use of Poisson–Nernst–Planck (PNP)
equations (Rubinstein 1990) remains an important tool for the investigation of
transport phenomena at scales below 100 nm (Schoch et al. 2008). Typically,
the systems of interest exhibit very small Reynolds numbers so that the Stokes
equation can be applied to describe the slow flow of the solvent such as water, for
example. However, there are very few geometries and systems where this highly
nonlinear set of governing equations can be solved analytically.
This paper presents a complete solution of the PNP–Stokes equations for a
specific system, namely, an infinite, cylindrical channel with negative, uniform
wall charges, containing an acidic solution of protons and water (figure 1). Co-ions
are not present.
This system is reminiscent of water-filled pores in polymer electrolyte
membranes (PEMs), as used in hydrogen fuel cells, which are proton exchange
membranes whose ion exchange mechanism between sulphonic acid groups and
water in nano-scale pores enables a current to flow through this charge-selective
electrolyte (Eikerling et al. 2001). Much research has been undertaken to better
*Author for correspondence (peter.berg@uoit.ca).

Received 28 January 2011

Accepted 17 May 2011 3157 This journal is © 2011 The Royal Society
 Downloaded from http://rspa.royalsocietypublishing.org/ on June 22, 2015

3158 P. Berg and J. Findlay

(a) – – – – – – – (b)
– –
r
R –
–
r
a
a

z – –
R

– –
–
– – – – – – –

Figure 1. Pore geometry.

understand the flow of water and protons in such PEM nano-pores (Zawodzinski
et al. 1995; Paddison et al. 1998, 2000a; Eikerling & Kornyshev 2001; Eikerling
et al. 2001).
Recently, SAXS spectra analysis has been used by Schmidt-Rohr & Chen
(2008) to postulate that the pore morphology of PEM consists mainly of
cylindrical nano-channels with a circular cross section. Based on this evidence,
Berg & Ladipo (2009) derived an analytical solution of Stokes flow in these
cylindrical pores, following the ideas of Qiao & Aluru (2003). The results go
beyond previous solutions of linearized Poisson–Boltzmann equations (Rice &
Whitehead 1965) but fail to include a pressure term in the Stokes equation, as
well as a migration term in the Nernst–Planck equation (Rubinstein 1990). In
other words, only advection is considered along the pore.
In a more comprehensive numerical study, Karimi & Li (2005) solved the
complete set of PNP–Stokes equations for a cylindrical PEM pore. Expressions
are derived for the proton and potential distributions, and the water flow. The
authors mention explicitly that the system cannot be solved analytically which is
true for a spatially varying permittivity as used in their work. However, we will
show in what follows that the entire system of equations can be solved in closed
form for a constant permittivity.
Therefore, this brief contribution will close an important gap and present a full
solution of the PNP equations, coupled with Stokes flow, for such channels. The
results also include expressions for water drag and conductivity, two quantities
of the system that can be compared with experimental data.
However, we would like to underline that our analysis does not include a
Stern layer that plays an important role at this scale. Consequently, the model
predictions might deviate from experimental results. Indeed, the focus of this
paper lies in the complete solution of the PNP–Stokes equations while the
reference to PEMs serves mainly as a motivation. It adds relevance to the results
but we do not claim that the analytical solution can fully describe PEM nano-
channels. To the contrary, there is strong evidence that a Stern layer must be
incorporated to reproduce experimental findings. For a discussion of the issues
surrounding the modelling of PEM nano-pores, their Stern layers and overlapping
electric double layers, we refer to Ladipo et al. (2011).

Proc. R. Soc. A (2011)

 Downloaded from http://rspa.royalsocietypublishing.org/ on June 22, 2015

Solution of PNP–Stokes equations 3159

The next section presents analytical expressions for the electric potential, the
proton concentration, the velocity profile, the mass flux, the current, the water
drag and the conductivity in cylindrical channels, described by PNP–Stokes
equations. The last section will conclude with a brief comparison of analytical
results to experimental data of PEMs.

2. Model and solution

Similar to the work of Berg & Ladipo (2009), we study an infinite, cylindrical and
circular pore with a uniform surface charge density, ss , for a given pore radius R
(figure 1). The flow is driven in the (axial) z-direction by an external electric
field1 and a pressure gradient. In case of a PEM nano-pore, we assume that all
protons are dissociated from the sulphonic acid groups.
Owing to the symmetry of the domain, the system equations can be split into
radial and axial parts. In the radial direction, the proton concentration, c, varies
as well as the electric potential, j, generated by them. The fluid velocity, ur , is
zero. In the axial direction, j and c do not vary, the fluid velocity, uz (r), (which
we will denote by u = u(r) in what follows) is non-zero and a constant external
electric field, E, is applied.
The system of equations includes the PNP equations for j and c
 
1 d dj qc
r =− , (2.1)
r dr dr ee0
VJ+ = 0 (2.2)
with  
dc(r) DF dj
J+ = qu(r)c(r)ê z − D ê r + qc(r) E ê z − ê r , (2.3)
dr RT dr
and a constant permittivity, e. In reality, the permittivity could be a function
of spatial coordinates or of system size (Paddison et al. 1998) but in order to
derive an analytical solution, it is a constant in this work. Also, owing to the
above assumptions and the symmetry of the flow, the axial contribution (ê z ) to
the divergence in equation (2.2) is zero.

(a) Electric potential and proton concentration

In the radial direction, J+ is zero and equations (2.1) and (2.3) can be combined
to yield the Poisson–Boltzmann equation (Berg & Ladipo 2009)
   
1 d dj qc0 qj
r =− exp − . (2.4)
r dr dr ee0 kB T
1 Note that we impose a very distinct flow scenario for this study. In principle, one does not need
to differentiate between internal and external fields. Rather, the field should be determined by
boundary conditions at the pore in- and outlet. Only particular geometries, e.g. a straight, infinite
and circular cylinder allow for the separation of the field into an ‘internal’ (radial) component and
an ‘external’ (axial) component.

Proc. R. Soc. A (2011)

 Downloaded from http://rspa.royalsocietypublishing.org/ on June 22, 2015

3160 P. Berg and J. Findlay

This nonlinear ordinary differential equation can be solved analytically for the
following boundary conditions (Berg & Ladipo 2009):

dj
=0 at r = 0 (2.5)
dr
and
dj ss
=− at r = R. (2.6)
dr ee0
The first equation is a symmetry condition. The last equation ensures global
electro-neutrality.2
Furthermore, we impose the gauge j(0) = 0 and the non-dimensional
coordinate3 x according to r = Rx. The derivation of the solutions is presented
by Berg & Ladipo (2009). They are
 
kB T l 2 2
j(x) = ln 1 − x (2.7)
q 8

and
c(x) = c0 (1 − lx 2 /8)−2 , (2.8)

where l is a dimensionless parameter

R 2 q 2 c0
l= , (2.9)
ee0 kB T
which relates to the Debye length

R ee0 kB T
d=√ = , (2.10)
l q 2 c0

with
8ss
c0 = . (2.11)
q 2 R2 s s /ee0 kB T − 4qR

Note that the solutions j(x) and c(x) do not include the diffusion coefficient
of protons since D scales out for the radial component of the Nernst–Planck
equation (2.3). In other words, both diffusion and mobility (i.e. migration) scale
proportional to D.

2 The expression (2.6) usually relates to the electric field of a charged plane. It is applicable in
our case since it can be derived for a small test volume at the wall whose curvature effectively
approaches zero in the limit of vanishing volume, thereby resembling a charged plane.
3 In this work, when we switch from r to x, the corresponding functions keep their symbolic

expressions for reasons of simplicity even though, strictly speaking, we should change notation.

Proc. R. Soc. A (2011)

 Downloaded from http://rspa.royalsocietypublishing.org/ on June 22, 2015

Solution of PNP–Stokes equations 3161

(b) Fluid velocity

The fluid velocity in the axial direction is given by the Stokes equation since
we are dealing with very slow flow (small Reynolds number, Re  1)
 
dp 1 d du
0=− +m r + qc(r)E. (2.12)
dz r dr dr
This can be solved analytically. First, we re-write this equation, employing the
non-dimensional coordinate x, as
 
m d du
x = −qc(x)ER2 + PR2 . (2.13)
x dx dx
Here, P is the constant pressure gradient dp/dz. The boundary conditions are
no-slip on the wall and a symmetry condition at the centre of the pore
u(1) = 0 (2.14)
and
du
=0 at x = 0. (2.15)
dx
In equation (2.13), substituting for c(x) by use of equation (2.8), and integrating
both sides leads to
 
du qc0 ER2 x PR2
x =− dx + x dx
dx m (1 − lx 2 /8)2 m
 
4qc0 ER2 lx PR2
=− dx + x dx (2.16)
ml 4(1 − lx 2 /8)2 m
du 4qc0 ER2 1 PR2 C1
⇒ =− + x+ . (2.17)
dx ml x(1 − lx /8)
2 2m x
The constant of integration is determined by the symmetry condition (2.15). In
order for the right-hand side to remain finite at x = 0, we need C1 = 4qc0 ER2 /lm,
leading to
 
du 4qc0 ER2 1 1 PR2
=− − + x
dx lm x(1 − lx 2 /8) x 2m
qc0 ER2 x PR2
=− + x. (2.18)
2m 1 − lx 2 /8 2m
Integrating both sides again, we find
u(x) = u0 (x) + up (x), (2.19)
 
2qc0 ER2 lx 2
with u0 (x) = ln 1 − + C2 (2.20)
lm 8
PR2 2
and up (x) = x + C3 , (2.21)
4m

Proc. R. Soc. A (2011)

 Downloaded from http://rspa.royalsocietypublishing.org/ on June 22, 2015

3162 P. Berg and J. Findlay

where u0 (x) is the solution found in Berg & Ladipo (2009) and up (x) is a solution
containing the pressure term. With the no-slip condition at the channel wall,
u(1) = 0, the integration constants and solutions are found to be
 
2qc0 ER2 l
C2 = − ln 1 −
lm 8
(2.22)
2
PR
and C3 = − .
4m
 ⎫
2qc0 ER2 1 − lx 2 /8 ⎪
⇒ u0 (x) = ln ,⎪
⎪
lm 1 − l/8 ⎬
(2.23)
PR2 ⎪
⎪
up (x) = − (1 − x ).
2 ⎪
⎭
4m
In summary, this results in a velocity profile
 
2qc0 ER2 1 − lx 2 /8 PR2
u(x) = ln − (1 − x 2 ). (2.24)
lm 1 − l/8 4m
We observe that the pressure term simply adds a Poiseuille-like component to
the electrically driven flow component.

(c) Total mass flux

To compute the water drag, we first need to determine the total mass flux, Ṁ ,
through a cross section of a cylinder. It is given by
r
Ṁ = 2p[ru(r)]r dr (2.25)
0
and 1
Ṁ = 2prR 2
x[u0 (x) + up (x)] dx. (2.26)
0
The first term is solved in Berg & Ladipo (2009) as follows:
1    
2qc0 ER2 1 4 l −1 1
xu0 (x) dx = ln 1 − − (2.27)
0 m l l 8 2
and hence 1
2qc0 ER2
xu0 (x) dx = f (l), (2.28)
0 m
where f (l) represents the expression in the brackets. The second term is found
to be
1 
PR2 1 PR2
xup (x) dx = − x(1 − x 2 ) dx = − . (2.29)
0 4m 0 16m
This gives the mass flux
 
2qc0 ER2 PR2
Ṁ = 2prR2 f (l) − (2.30)
m 16m

Proc. R. Soc. A (2011)

 Downloaded from http://rspa.royalsocietypublishing.org/ on June 22, 2015

Solution of PNP–Stokes equations 3163

and therefore
4prqc0 ER4 prPR4
Ṁ = f (l) − , (2.31)
m 8m
which scales like R4 , similar to Poiseuille flow.

(d) Conductivity and water drag

Following the approach in Berg & Ladipo (2009), the model predictions can be
compared with experimental results for conductivity and water drag in PEMs.
To do so, the total proton flux, or current, is required which is determined by
1  
DF
I = 2pR qc(x)u(x) +
2
qEc(x) x dx (2.32)
0 RT
or
1  
DF
I = 2pR2 qxc(x)[u0 (x) + up (x)] + qExc(x) dx. (2.33)
0 RT
The integration will be split into three steps. The first step is completed in Berg &
Ladipo (2009), involving the u0 (x) term
1 
2qc02 ER2 1 ln(1 − lx 2 /8) − ln(1 − l/8)
c(x)u0 (x)x dx = dx
0 lm 0 (1 − lx 2 /8)2
  
2qc02 ER2 1 4 l
= + ln 1 − (2.34)
m 2l(1 − l/8) l2 8
or
1
2qc02 ER2
c(x)u0 (x)x dx = g(l), (2.35)
0 m
where g(l) represents the term in the brackets. The second step involves the
solution up (x), which includes the pressure
1 1  
c0 PR2
c(x)up (x)x dx = − (1 − x ) dx
2
0 (1 − lx /8)
2 2 4m
0

c0 PR2 1 x
=− (1 − x 2 ) dx (2.36)
0 (1 − lx /8)
4m 2 2

and therefore
1  1 
c0 PR2 1 x x
c(x)up (x)x dx = − dx − x 2
dx , (2.37)
0 (1 − lx /8) 0 (1 − lx /8)
4m 2 2 2 2
0

where each integral is solved separately. The first integral is similar to an integral
found in equation (2.16), namely
1  
x 4 1
dx = − 1 . (2.38)
0 (1 − lx /8) l 1 − l/8
2 2

Proc. R. Soc. A (2011)

 Downloaded from http://rspa.royalsocietypublishing.org/ on June 22, 2015

3164 P. Berg and J. Findlay

The second term can be solved using integration by parts as follows:

1   1  1
x 2 4 1 4 1
x 2
dx = x − 2x dx (2.39)
0 (1 − lx 2 /8)2 l 1 − lx 2 /8
0 0 l (1 − lx 2 /8)

which yields
1    
x 4 1 32 l
x 2 dx = + 2 ln 1 − . (2.40)
0 (1 − lx /8)
2 2 l 1 − l/8 l 8
Combining equations (2.38) and (2.40) gives
1       
c0 PR2 4 1 4 1 32 l
c(x)up (x)x dx = − −1 − − 2 ln 1 −
0 4m l 1 − l/8 l 1 − l/8 l 8
  
c0 PR2 4 32 l
=− − − 2 ln 1 − (2.41)
4m l l 8
or
1
c0 PR2
c(x)up (x)x dx = − h(l), (2.42)
0 4m
where h(l) represents the term in the brackets. The final term in equation (2.33)
is found to be
 
2pR2 DFqE 1 2pR2 DFqc0 E 1 x
c(x)x dx = dx
0 (1 − lx /8)
RT RT 2 2
0
  
2pR2 DFqc0 E 4 1
= −1 (2.43)
RT l 1 − l/8
and hence
1
2pR2 DFqE 2pDFqc0 E
c(x)x dx = k(l)R2 , (2.44)
RT 0 RT
where k(l) represents the term in the brackets. Finally, combining
equations (2.35), (2.42) and (2.44) gives the solution for the total proton flux

4pq 2 c02 E pqc0 P 2pDFqc0 E

I= g(l)R4 − h(l)R4 + k(l)R2 . (2.45)
m 2m RT

Note that in contrast to Ṁ , I does not exhibit a R4 -scaling.

The conductivity, sp , can now be computed as (Berg & Ladipo 2009)

I 4q 2 c02 qc0 P 2DFqc0

sp = = g(l)R2 − h(l)R2 + k(l). (2.46)
pR2 E m 2mE RT
It is an expression that depends trivially on all system parameters except for the
Debye length. However, it does not scale like R2 . In fact, it contains a term (the
last term, representing migration) that does not scale with R at all.

Proc. R. Soc. A (2011)

 Downloaded from http://rspa.royalsocietypublishing.org/ on June 22, 2015

Solution of PNP–Stokes equations 3165

The second quantity to be computed is the electro-osmotic drag, hdrag . As in

Berg & Ladipo (2009), it is given by the ratio of the water flux to proton flux.
Both quantities are converted to molar fluxes by including Faraday’s constant
and the molar weight of water, w = 0.018 kg mol−1 , to yield
Ṁ F
hdrag (R) =
· . (2.47)
w I
It depends trivially on all system parameters except for the Debye length through
the functions f , g, h and k. Again, there is no scaling with R.
In fact, neither sp nor hdrag exhibit a Ra -scaling, for some real number a.

3. Results and discussion

The analytical results above are very helpful to obtain order of magnitude
estimates for water and proton flows in (PEM) nano-channels. In addition, they
can play an important role in validating numerical models of such systems for
more complex three-dimensional geometries. Computer codes can be validated
first for straight-channel geometries before they are modified to simulate more
intricate domains. For this reason alone, one could argue that the above results
are an important contribution to the field of electro-kinetic theory.

(a) Proton and water flow in polymer electrolyte membranes

Notwithstanding the limited applicability of the PNP–Stokes equations to
PEM nano-channels, especially in light of neglecting a Stern layer and a varying
permittivity, we will now briefly discuss two characteristic quantities of PEM that
can be measured experimentally.
Figures 2 and 3 show the dependency of the water drag and the conductivity
on the radius of the pore for different values of the diffusion coefficient and the
pressure gradient. The latter opposes the motion of protons just like in a PEM
fuel cell. Here, back diffusion can occur when the cathode is sufficiently hydrated,
compared with the anode. The required system parameters, listed in table 1, are
adopted from the study by Berg & Ladipo (2009).
Water drag tends to increase with R as more water flows through the pore
while the total amount of protons remains fixed (since ss is assumed to scale like
1/R, thereby simulating simple pore swelling and conservation of wall charges).
However, at large R, hdrag decreases eventually since the pressure starts to have
a larger effect on water flow as R grows, which is a feature reminiscent of
Poiseuille flow.
We see that only for very large pressure gradients of about 10 × 1010 Pa m−1 ,
corresponding to about 100 atmospheres across a 100 mm PEM, the pressure has
an impact on the water drag. Ultimately, the drag coefficient can turn negative,
equivalent to a water flow opposing the proton flow. This is referred to as back
diffusion in the PEM fuel cell literature. Such large pressure gradients, required to
observe back diffusion, would have to be generated in a PEM fuel cell by capillary
effects at the membrane surface, specifically at the pore necks, when the cathode
is fully hydrated and the anode is under-saturated (Eikerling & Berg 2011). This
will be subject to future investigation but it holds the promise of explaining back
diffusion quantitatively.

Proc. R. Soc. A (2011)

 Downloaded from http://rspa.royalsocietypublishing.org/ on June 22, 2015

3166 P. Berg and J. Findlay

D = 0, P = 0
D = 7.5 × 10−10, P = 0
D = 7.5 × 10−10, P = 1.0 × 1010
D = 7.5 × 10−10, P = 3.0 × 1010
D = 7.5 × 10−10, P = 5.0 × 1010
D = 7.5 × 10−10, P = 10.0 × 1010

60

50

40

30
hdrag

20

10

−10
0 0.5 1.0 1.5 2.0 2.5
R (× 10–9)

Figure 2. Water drag as a function of pore radius R (in nanometres) for various values of the
diffusion coefficient (in m2 s−1 ) and the pressure gradient (in Pa m−1 ).

80
D = 0, P = 0
D = 7.5 × 10−10, P = 0
60 D = 7.5 × 10−10, P = 1.0 × 1010
D = 7.5 × 10−10, P = 3.0 × 1010
sp (S m–1)

D = 7.5 × 10−10, P = 5.0 × 1010

40
D = 7.5 × 10−10, P = 10.0 × 1010

20

0 0.5 1.0 1.5 2.0 2.5

R (× 10–9)

Figure 3. Pore conductivity as a function of pore radius R (in nanometres) for various values of
the diffusion coefficient (in m2 s−1 ) and the pressure gradient (in Pa m−1 ).

Meanwhile, the water drag coefficient at R = 1 nm is roughly 10 when migration

and a pressure gradient are included and, therefore, still too large when compared
with experimental data (Zawodzinski et al. 1995) which indicates hdrag ≈ 1 − 3.
Strictly speaking, hdrag should be measured when P = 0 (the two top curves),

Proc. R. Soc. A (2011)

 Downloaded from http://rspa.royalsocietypublishing.org/ on June 22, 2015

Solution of PNP–Stokes equations 3167

Table 1. Parameters and variables of the system.

variable,
parameter description value

c0 proton concentration at pore centre

D diffusion coefficient of protons
E applied external electric field 0.1 V/(50 × 10−6 m) = 2000 V m−1
F Faraday constant
kB Boltzmann constant
P pressure gradient
q electron charge
R radius
R ideal gas constant
T temperature 353 K
u fluid velocity
w molar weight of water
e permittivity 45
e0 vacuum permittivity
m viscosity 3.35 × 10−4 Pa s
r density of water
ss surface charge density −q/(10−9 m)2 = 0.16 C m−2

making the result look even worse. In this sense, our simple continuum model
for the electro-kinetic flow in a PEM pore fails. However, the results might
change substantially when a non-constant permittivity is included, as observed
in microwave experiments of Nafion (Paddison et al. 2000b) and theoretical
calculations (Paul & Paddison 2004a,b), along with a larger viscosity typical for
PEM nano-pores, a Stern layer and a pore network approach (see Ladipo et al.
2011, §2f ). Also, the value for the proton diffusion coefficient, D, remains a topic
of controversy in the literature.
In simulations of a PNP–Stokes system for a Nafion pore, Ladipo et al. (2011)
have shown that a non-constant permittivity, increasing from about 10 near the
pore wall to the value of bulk water (80) near the pore centre (Paul & Paddison
2004a,b), can indeed reduce the water drag by a factor of 2. Yet, only the inclusion
of a Stern layer, containing immobile water molecules, addresses the water drag
discrepancy successfully, bringing its value in line with experimental data.
On the other side, the conductivity exhibits reasonable values that are on the
right order of magnitude when compared with fully saturated Nafion membranes
(sp = 0.07 S cm−1 ). With a vanishing pressure gradient, we find a value of
sp ≈ 0.20 S cm−1 , about 50 per cent larger than the reference value of Berg &
Ladipo (2009), i.e. sp ≈ 0.12 S cm−1 , for a radius R = 1 nm. Pressure (via P) only
has a significant impact on the results at extremely large gradients. Further
modifications to the computation of an effective conductivity would have to
include a re-scaling that incorporates the size of the backbone and, again, an
upscaling through a pore network approach. The inclusion of the backbone alone
reduces the conductivity by a factor of two (sp ≈ 0.10 S cm−1 ), thereby predicting
a very reasonable value close to experimental data.

Proc. R. Soc. A (2011)

 Downloaded from http://rspa.royalsocietypublishing.org/ on June 22, 2015

3168 P. Berg and J. Findlay

Finally, it should be noted that discrepancies between theoretical and

experimental results may be explained by the inherent limitations of the PNP–
Stokes modelling approach that is usually applied to dilute solutions. In contrast,
a PEM nano-pore contains a very high concentration of protons. These, as well
as the water molecules, interact strongly with both the negative side chains
(Paddison et al. 2000b) and the polymer backbone (i.e. the pore walls). It is
a delicate matter to include such intricate details into continuum models. One
alternative might be the employment of Stefan–Maxwell equations in which the
interaction between water, protons and side chains can be incorporated in an
explicit manner. This allows for some level of flexibility when modelling water
and proton flows, unlike the force balance in the Stokes equation (2.12). Here,
the force that protons exert onto water, qc(r)E, equals the force that the electric
field exerts onto the protons. In other words, this rigid interaction between water
and protons does not reflect the rich phenomena of proton transport mechanisms
in these pores, especially as a function of the water content and the distance from
the side chains (Paul & Paddison 2004b).
The authors would like to thank Toyota Motor Engineering and Manufacturing North America
(TEMA) for financial support of this research.

References
Berg, P. & Ladipo, K. 2009 Exact solution of an electro-osmotic flow problem in a cylindrical
channel of polymer electrolyte membranes. Proc. R. Soc. A 465, 2663–2679. (doi:10.1098/rspa.
2009.0067)
Eikerling, M. & Berg, P. 2011 Poroelectroelastic theory of water sorption and swelling in polymer
electrolyte membranes. Soft Matter 1–16. (doi:10.1039/c1sm05273j)
Eikerling, M. & Kornyshev, A. 2001 Proton transfer in a single pore of a polymer electrolyte
membrane. J. Electroanal. Chem. 502, 1–14. (doi:10.1016/S0022-0728(00)00368-5)
Eikerling, M., Kornyshev, A., Kuznetsov, A., Ulstrup, J. & Walbran, S. 2001 Mechanisms of proton
conductance in polymer electrolyte membranes. J. Phys. Chem. B 105, 3646–3662. (doi:10.1021/
jp003182s)
Karimi, G. & Li, X. 2005 Electroosmotic flow through polymer electrolyte membranes in PEM fuel
cells. J. Power Sources 140, 1–11. (doi:10.1016/j.jpowsour.2004.08.018)
Ladipo, K. O., Berg, P., Kimmerle, S. J. & Novruzi, A. 2011 Effects of radially dependent
parameters on proton transport in PEM nanopores. J. Chem. Phys. 134, 074103. (doi:10.1063/
1.3552232)
Paddison, S., Reagor, D. & Zawodzinski, T. 1998 High frequency dielectric studies of hydrated
Nafion. J. Electroanal. Chem. 459, 91–97. (doi:10.1016/S0022-0728(98)00321-0)
Paddison, S., Paul, R. & Zawodzinski, T. 2000a A statistical mechanical model of proton and water
transport in a proton exchange membrane. J. Electrochem. Soc. 147, 617–626. (doi:10.1149/
1.1393243)
Paddison, S., Bender, G., Kreuer, K. D., Nicoloso, N. & Zawodzinski, T. A. 2000b The microwave
region of the dielectric spectrum of hydrated Nafion (R) and other sulfonated membranes. J. New
Mat. Electrochem. Sys. 3, 291–300.
Paul, R. & Paddison, S. 2004a The phenomena of dielectric saturation in the water domains of
polymer electrolyte membranes. Solid State Ionics 168, 245–248. (doi:10.1016/j.ssi.2003.06.001)
Paul, R. & Paddison, S. 2004b Structure and dielectric saturation of water in hydrated polymer
electrolyte membranes: inclusion of the internal field energy. J. Phys. Chem. B 108, 13 231–
13 241. (doi:10.1021/jp048501k)
Qiao, R. & Aluru, N. 2003 Ion concentrations and velocity profiles in nanochannel electroosmotic
flows. J. Chem. Phys. 118, 4692–4701. (doi:10.1063/1.1543140)

Proc. R. Soc. A (2011)

 Downloaded from http://rspa.royalsocietypublishing.org/ on June 22, 2015

Solution of PNP–Stokes equations 3169

Rice, C. & Whitehead, R. 1965 Electrokinetic flow in a narrow cylindrical capillary. J. Phys. Chem.
69, 4017–4024. (doi:10.1021/j100895a062)
Rubinstein, I. 1990 Electro-diffusion of ions, vol. I, pp. 23–56. Philadelphia, PA: SIAM.
Schmidt-Rohr, K. & Chen, Q. 2008 Parallel cylindrical water nanochannels in Nafion fuel-cell
membranes. Nat. Mater. 7, 75–83. (doi:10.1038/nmat2074)
Schoch, R., Han, J. & Renaud, P. 2008 Transport phenomena in nanofluidics. Rev. Mod. Phys. 80,
839–883. (doi:10.1103/RevModPhys.80.839)
Zawodzinski, T., Davey, J., Valerio, J. & Gottesfeld, S. 1995 The water content dependence
of electro-osmotic drag in proton-conducting polymer electrolytes. Electrochem. Acta
40, 297–302. (doi:10.1016/0013-4686(94)00277-8)

Proc. R. Soc. A (2011)