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.
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 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).
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, σs, 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, ψ, generated by them. The fluid velocity, ur, is zero. In the axial direction, ψ 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 ψ and c 2.1 2.2 with 2.3 and a constant permittivity, ϵ. 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 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) 2.4 This nonlinear ordinary differential equation can be solved analytically for the following boundary conditions (Berg & Ladipo 2009): 2.5 and 2.6 The first equation is a symmetry condition. The last equation ensures global electro-neutrality.2
Furthermore, we impose the gauge ψ(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 2.7 and 2.8 where λ is a dimensionless parameter 2.9 which relates to the Debye length 2.10 with 2.11 Note that the solutions ψ(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.
(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) 2.12 This can be solved analytically. First, we re-write this equation, employing the non-dimensional coordinate x, as 2.13 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 2.14 and 2.15 In equation (2.13), substituting for c(x) by use of equation (2.8), and integrating both sides leads to 2.16 2.17 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=4qc0ER2/λμ, leading to 2.18 Integrating both sides again, we find 2.19 2.20 and 2.21 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 2.22 2.23 In summary, this results in a velocity profile 2.24 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 2.25 and 2.26 The first term is solved in Berg & Ladipo (2009) as follows: 2.27 and hence 2.28 where f(λ) represents the expression in the brackets. The second term is found to be 2.29 This gives the mass flux 2.30 and therefore 2.31 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 2.32 or 2.33 The integration will be split into three steps. The first step is completed in Berg & Ladipo (2009), involving the u0(x) term 2.34 or 2.35 where g(λ) represents the term in the brackets. The second step involves the solution up(x), which includes the pressure 2.36 and therefore 2.37 where each integral is solved separately. The first integral is similar to an integral found in equation (2.16), namely 2.38 The second term can be solved using integration by parts as follows: 2.39 which yields 2.40 Combining equations (2.38) and (2.40) gives 2.41 or 2.42 where h(λ) represents the term in the brackets. The final term in equation (2.33) is found to be 2.43 and hence 2.44 where k(λ) represents the term in the brackets. Finally, combining equations (2.35), (2.42) and (2.44) gives the solution for the total proton flux 2.45 Note that in contrast to , I does not exhibit a R4-scaling.
The conductivity, σp, can now be computed as (Berg & Ladipo 2009) 2.46 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.
The second quantity to be computed is the electro-osmotic drag, ηdrag. 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 2.47 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 σp nor ηdrag exhibit a Rα-scaling, for some real number α.
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 σs is assumed to scale like 1/R, thereby simulating simple pore swelling and conservation of wall charges). However, at large R, ηdrag 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 μm 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.
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 ηdrag≈1−3. Strictly speaking, ηdrag should be measured when P=0 (the two top curves), 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 (σp=0.07 S cm−1). With a vanishing pressure gradient, we find a value of σp≈0.20 S cm−1, about 50 per cent larger than the reference value of Berg & Ladipo (2009), i.e. σp≈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 (σp≈0.10 S cm−1), thereby predicting a very reasonable value close to experimental data.
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.
↵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.
↵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.
- Received January 28, 2011.
- Accepted May 17, 2011.
- This journal is © 2011 The Royal Society