Induced charge electroosmotic flow with finite ion size and solvation effects

TL;DR

This paper introduces a new numerical method for modeling electroosmotic flows that accounts for finite ion size and solvation effects, improving accuracy in boundary layer simulations of electrolyte systems.

Contribution

It develops a thermodynamically consistent numerical approach that incorporates volume exclusion and solvation effects into electroosmotic flow modeling.

Findings

Successfully simulates vortex generation in nanochannels due to electroosmotic forces.

Ensures numerical stability and thermodynamic consistency in boundary layer modeling.

Demonstrates the importance of finite ion size effects in electroosmotic flow simulations.

Abstract

The most common mathematical models for electrolyte flows are based on the dilute solution assumption, leading to a coupled system of the Nernst--Planck--Poisson drift-diffusion equations for ion transport and the Stokes resp. Navier--Stokes equations for fluid flow of the solvent. In charged boundary layers the dilute solution assumption is in general not valid and volume exclusion and solvation effects have to be taken into account in a thermodynamically consistent way. Whenever boundary layer effects have a dominant impact on the global behavior of a certain electrochemical system, an accurate numerical simulation depends on the correct incorporation of these effects. In this contribution we present a novel numerical solution approach which aims at preserving on the discrete level consistency with basic thermodynamic principles and structural properties like independence of flow velocities from gradient contributions to external forces. We illustrate capabilities of the method by an example of vortex generation due to induced charge electroosmotic forces at an electrode inside a nanochannel.

Figures

Equations

$$\rho {\partial }_t\vec{u}-\nu \Delta \vec{u}+\rho (\vec{u}\cdot \nabla )\vec{u}+\nabla p$$

$$\nabla \cdot (\rho \vec{u})$$

$${\partial }_t{c}_i+\nabla \cdot ({\vec{N}}_i+c_i\vec{u})$$

$$-\nabla \cdot (\epsilon \nabla \psi )$$

$${\vec{N}}_i$$

$${\mu }_i$$

$$1$$

$$\rho$$

$${\vec{N}}_i$$

$${\mu }_i$$

$$1$$

$$\rho$$

$${\vec{N}}_i$$

$$-\nabla \cdot (\epsilon \nabla \psi )$$

$$0$$

$${\vec{N}}_i$$

$$a_i$$

$${\beta }_i=\beta =$$

$$0\le c_i={\beta }_i{a}_i\le \frac{1}{v_0}.$$

$${\beta }_i$$

$${\partial }_{t}c+\nabla \cdot \vec{J}=0$$

$$0=$$

$$0=|K|\left(c_K^n-c_K^{n-1}\right)+\sum _{\begin{array}{c}L\text{neighbor of}K\end{array}}\left(t^n-t^{n-1}\right)|\partial K\cap \partial L|J_{KL}^n,$$

$$J_{KL}=D{\beta }_{KL}\cdot \frac{B\left(-{\delta }_{KL}-\frac{u_{KL}}{D}\right)a_K-B\left({\delta }_{KL}+\frac{u_{KL}}{D}\right)a_L}{h_{KL}},$$

$$u_{KL}$$

$$\sum _{\begin{array}{c}L\text{neighbor of}K\end{array}}|\partial K\cap \partial L|u_{KL}=0$$

$${\int }_{\Omega }\nu \nabla {\vec{u}}_h:\nabla {\vec{v}}_{h}dx+{\int }_{\Omega }{p}_h\nabla \cdot {\vec{v}}_{h}dx$$

$${\int }_{\Omega }{q}_h\nabla \cdot {\vec{u}}_{h}dx$$

$${\int }_{\Omega }{q}_h\nabla \cdot {\vec{v}}_{h}dx=0\text{for all }{q}_h\in Q_h⟹\nabla \cdot (\Pi {\vec{v}}_h)=0\text{in }\Omega .$$

Full text

1 Introduction

With the miniaturization of electrochemical devices and the trend towards nanometer scale electrochemistry [1, 2, 3], understanding the behavior of liquid electrolytes at the micro- and nanoscale becomes increasingly important. Of particular interest are correct double layer models, the behavior of electroosmotic flows at the nanoscale and their interaction with species transport. Due to the complex physical interactions present in electrolyte flows, numerical simulation techniques are important tools for developing a deeper understanding of the flow behavior.

This contribution introduces a simulation approach for coupled ion transport and electrolyte flow. Section 2 introduces a modified Nernst–Planck–Poisson–Navier–Stokes model which has its foundations in first principles of nonequilibrium thermodynamics [4, 5, 6] and takes into account ion-solvent interactions, finite ion size and solvation effects.

Section 3 introduces a finite volume discretization approach for ion transport in a self-consistent electric field which is motivated by results from semiconductor device simulation [7]. Pressure robust mixed finite element methods for fluid flow [8], and a fix point approach for coupling to ion transport are introduced.

Section 4 provides a numerical example which shows the potential significance of the proposed approach for modeling electrochemical processes at the nanoscale.

2 The model

Electroosmotic flows are characterized by the presence of an electric field that exerts a net force on the fluid molecules in regions where the local net charge due to the present ions is nonzero. Being part of the momentum balance for the barycentric velocity of the fluid, this net force induces fluid motion. Correspondingly, ions are advected by the barycentric velocity field. Their motion relative to the barycentric velocity is induced by the gradients of their chemical potential and the electrostatic potential. The elastic interactions between the ions and the solvent exerts a counterforce against ion motion.

The ion charge density varies with the redistribution of the ions and contributes to the electric field.

The dilute solution assumption lying at the foundation of classical electrolyte models [9, 10] assumes that the ion concentration in the electrolyte is small compared to the concentration of the solvent and consequently neglects the ion-solvent interaction. However, accumulation of ions in polarization boundary layers violates this assumption. As a consequence, the resulting electrolyte model e.g. severely overestimates double layer capacitances at ideally polarizable electrodes [11].

The model used throughout this paper has been introduced in [4, 5, 6] based on consistent application of the principles of nonequilibrium thermodynamics [12]. It includes ion volume and solvation effects and consistently couples the transport equations to the momentum balance. It generalizes previous approaches to include steric (ion size) effects [13, 14, 15, 16, 17], see also [18].

2.1 Model equations

In a given bounded space-time domain $\Omega\times(0,t^{\sharp})\subset\mathbb{R}^{d}\times(0,\infty)$, and with appropriate initial and boundary conditions, system (1)–(4) describes the isothermal evolution of the molar concentration of $N$ charged species with molar densities (concentrations) $c_{1}\dots c_{N}$ with charge numbers $z_{1}\dots z_{N}$ dissolved in a solvent of concentration $c_{0}$. Species are further characterized by their molar volumes $v_{i}$ and molar masses $M_{i}$. The electric field is described as the gradient of the electrostatic potential $\psi$. The barycentric velocity of the mixture is denoted by $\vec{u}$, and $p$ is the pressure. The following equations are considered:

[TABLE]

Equation (1a) together with (1b) comprises the incompressible Navier–Stokes equations for a fluid of viscosity $\nu$ and constant density $\rho$. In the general case, where molar volumes and molar masses are not equal, $\rho$ would depend on the local composition of the electrolyte.

In regions where the charge density $q=F\sum_{i=1}^{N}z_{i}c_{i}$ ($F$ being the Faraday constant) is nonzero, the electric field $-\nabla\psi$ becomes a driving force of the flow.

The partial mass balance equations (1c) describe the redistribution of species concentrations due to advection in the velocity field $\vec{u}$ and molar diffusion fluxes $\vec{N}_{i}$. The Poisson equation (1d) describes the distribution of the electrostatic potential $\psi$ under a given configuration of the charge density $q$. The constant $\varepsilon$ is the dielectric permittivity of the medium.

The fluxes $\vec{N}_{i}$, the molar chemical potentials $\mu_{i}$ and the incompressibility constraint for a liquid electrolyte are given by

[TABLE]

The modified Nernst–Planck flux (2a) combines the gradients of the species chemical potentials $\nabla\mu_{i}$, the gradient of the solvent chemical potential $\nabla\mu_{0}$ and the electric field $-\nabla\psi$ as driving forces for the motion of ions of species $i$ relative to the barycentric velocity $\vec{u}$. In this equation, $D_{i}$ are the diffusion coefficients, $R$ is the molar gas constant, and $T$ is the temperature. Equation (2b) is a constitutive relation for the chemical potential $\mu_{i}$ depending on the local pressure and concentration. Here, $p^{\circ}$ is a reference pressure, and $\overline{c}=\sum_{i=0}^{N}c_{i}$ is the total species concentration. In (2c) a simple model for solvated ions is applied, see [5, 7, 19], which describes the mass and volume of a solvated ion by $\kappa_{i}M_{0}+M_{i}$ and $\kappa_{i}v_{0}+v_{i}$, respectively. The incompressibility constraint (2c) limits the accumulation of ions in the polarization boundary layer to physically reasonable values [5, 7].

The mass density of the mixture is

[TABLE]

As for reasonable solvation numbers $\kappa_{i}\approx 10$, molar masses and molar volumes of the solvated ions are dominated by the solvent mass and volume. Therefore we assume, for simplicity, $(\kappa_{i}M_{0}+M_{i})\approx(\kappa_{i}+1)M_{0}$ and $(\kappa_{i}v_{0}+v_{i})\approx(\kappa_{i}+1)v_{0}$. This approximation leads to a simplified constitutive function and, in particular, to a constant mass density of the mixture,

[TABLE]

Comparing the constitutive equations (4a)-(4c) to the classical Nernst–Planck flux [9, 10]

[TABLE]

which considers dilute solutions, one observes that in (5) the ion-solvent interaction described by the term $\nabla\mu_{0}$ is ignored. Moreover in (5) implicitly a material model is assumed that neglects the pressure dependence of $\mu_{i}$, which is inappropriate in charged boundary layers [4].

2.2 Reformulation in species activities

In order to develop a space discretization approach for the Nernst–Planck fluxes (4a), after [20], the system is reformulated in terms of (effective) species activities $a_{i}=\exp\left(\frac{\mu_{i}-(\kappa_{i}+1)\mu_{0}}{RT}\right)$. The quantity $\mu_{i}-(\kappa_{i}+1)\mu_{0}$ is sometimes denoted as entropy variable [21]. Introducing the activity coefficients $\gamma_{i}=\frac{a_{i}}{c_{i}}$ and its inverse (reciprocal) $\beta_{i}=\frac{1}{\gamma_{i}}=\frac{c_{i}}{a_{i}}$ allows to transform the Nernst–Planck–Poisson system consisting of (1c), (1d), (4a) to

[TABLE]

From (4b) and (4c) one obtains

[TABLE]

which is a linear system of equations which allows to express $\beta_{1}\dots\beta_{n}$ through $a_{1}\dots a_{n}$. It has the unique solution [20]

[TABLE]

It follows immediately that for any nonnegative solution $a_{1}\dots a_{n}$ of system (6), the resulting concentrations are bounded in a physically meaningful way:

[TABLE]

In the general case with different molar volumes and molar masses, system (7) becomes nonlinear, the quantities $\beta_{i}$ differ between species and in addition depend on the pressure $p$ [20, 7], leading to a nonlinear system of equations

[TABLE]

3 A coupled Finite-Volume-Finite-Element discretization

3.1 Thermodynamically consistent finite volume methods for species transport

A two point flux finite volume method on boundary conforming Delaunay meshes is used to approximate the Nernst–Planck–Poisson part of the problem. It has been inspired by the successful Scharfetter–Gummel box method for the solution of charge transport problems in semiconductors [22, 23]. For a recent overview on this method see [24]. It was initially developed for drift-diffusion problems in non-degenerate semiconductors exhibiting Boltzmann statistics for charge carrier densities whose fluxes are equivalent to the classical Nernst–Planck flux (5).

The time axis is subdivided into intervals $0=t^{0}<t^{1}<\dots<t^{N_{t}}=t^{\sharp}.$The simulation domain $\Omega$ is partitioned into a finite number of closed convex polyhedral control volumes $K\in\mathcal{K}$ such that $K\cap L=\partial K\cap\partial L$ and $\overline{\Omega}=\bigcup_{K\in\mathcal{K}}K$. With each control volume a node $\vec{x}_{K}\in K$ is associated. If the control volume intersects with the boundary $\partial\Omega$, its corresponding node shall be situated on the boundary: $\vec{x}_{K}\in\partial\Omega\cap K$. The partition shall be admissible [25], that is for two neighboring control volumes $K,L$, the edge $\overline{\vec{x}_{K}\vec{x}_{L}}$ is orthogonal to the interface between the control volumes $\partial K\cap\partial L$. Let $\vec{h}_{KL}=\vec{x}_{L}-\vec{x}_{K}$ and $h_{KL}=|\vec{h}_{KL}|$. Then, the normal vectors to $\partial K$ can be calculated as $\vec{n}_{KL}=\frac{1}{h_{KL}}\vec{h}_{KL}$.

A constructive way to obtain such a partition is based on the creation of a boundary conforming Delaunay triangulation resp. tetrahedralization of the domain and the subsequent construction of its dual, the Voronoi tessellation intersected with the domain, see e.g. [23, 26, 24], see also Figure 1.

Denote by $\vec{J}_{i}=c_{i}\vec{u}+\vec{N}_{i}=\beta_{i}a_{i}\vec{u}+\vec{N}_{i}$ the convection diffusion flux of the model under consideration. The general approach to derive a two point flux finite volume scheme for a conservation law (index $i$ omitted)

[TABLE]

consists in integrating the equation over a space-time control volume $K\times[t^{n-1},t^{n}]$:

[TABLE]

This is approximated via

[TABLE]

and it remains to define the numerical fluxes $J^{n}_{KL}$ which should approximate the continuous fluxes between two neighboring control volumes and depend on the unknown values in the two collocation points $\vec{x}_{K}$ and $\vec{x}_{L}$ at moment $t^{n}$ in order to obtain a fully implicit in time scheme.

The modification of the Scharfetter-Gummel scheme [22] proposed in [20] is based on the similarity of the expressions (6c) and (5). The flux expression (5) is the same as the drift-diffusion flux in non-degenerate semiconductors, for which this discretization scheme was initially derived. The only difference between (6c) and (5) is the multiplication with the pre-factor $\beta$. Therefore it appears to be reasonable to mimic this structure at the discrete level, and to derive a discrete equivalent of (6c) from the discrete version of (5) by multiplication with a suitable average of $\beta$: set

[TABLE]

where $B(\xi)=\frac{\xi}{e^{\xi}-1}$ be the Bernoulli function. and $\beta_{KL}$ is an average of the inverse activity coefficients $\beta_{K}$ and $\beta_{L}$, $\delta_{KL}=\frac{zF}{RT}(\psi_{K}-\psi_{L})$ is proportional to the local electric force, and

[TABLE]

is the normal integral over the interface $\partial K\cap\partial L$ of the convective flux scaled by $h_{KL}$. If the continuous flux is divergence free, i.e. it fulfills equation (1b), the flux projections $u_{KL}$ fulfill the discrete divergence condition

[TABLE]

which in the absence of charges and coupling through the activity coefficients guarantees a discrete maximum principle of the approximate solution [27].

The resulting time discrete implicit Euler finite volume upwind scheme guarantees nonnegativity of discrete activities and exact zero fluxes under thermodynamic equilibrium conditions. Moreover it guarantees the bounds (8) [20].

3.2 Pressure robust, divergence free finite elements for fluid flow.

Mixed finite element methods approximate the Stokes resp. Navier–Stokes equations based on an inf-sup stable pair of velocity ansatz space $V_{h}$ and pressure ansatz space $Q_{h}$. A fundamental property of the Stokes and Navier–Stokes equations consists in the fact that — under appropriate boundary conditions — the addition of a gradient force to the body force on the right-hand side of the momentum balance (1a) leaves the velocity unchanged, as it just can be compensated by a change in the pressure. Most classical mixed finite element methods for the Navier–Stokes equations do not preserve this property [28]. As a consequence, the corresponding error estimates for the velocity depend on the pressure [29]. Moreover, the discrete solution $\vec{u}_{h}$ fulfills the divergence constraint only in a discrete finite element sense. This raises problems when coupling the flow simulation to a transport simulation using finite volume methods, because the maximum principle for the species concentration is directly linked to the divergence constraint in the finite volume sense (11) [27] that is different to the finite element sense.

The pressure robust mixed methods first introduced in [8] and comprehensively described in [30], are based on the introduction of a divergence free velocity reconstruction operator $\Pi$ into the discrete weak formulation of the flow problem. The resulting discretization of the stationary Stokes equation (provided here for simplicity) reads as: find $(\vec{u}_{h},p_{h})\in V_{h}\times Q_{h}$ such that

[TABLE]

This formulation differs from that of the classical mixed methods only in the right hand side. The reconstruction operator $\Pi$ has the following properties:

1. (i)

If $\vec{v}_{h}\in V_{h}$ is divergence free in the weak sense, then its reconstruction $\Pi\vec{v}_{h}$ is pointwise divergence free:

[TABLE] 2. (ii)

The change of the test function by the reconstruction operator causes a consistency error that should have the same asymptotic convergence rate of the original method and should not depend on the pressure.

Under these conditions, the resulting velocity error estimate is independent of the pressure [30] and has the optimal expected convergence rate. Furthermore, a good velocity approximation can be obtained without the need to resort to high order pressure approximations. This allows significant reduction of degrees of freedom that are necessary to obtain a prescribed accuracy of the velocity. The action of $\Pi$ on a discrete velocity field can be calculated locally, on elements or element patches. Therefore its implementation leads to low overhead in calculations.

3.3 Coupling strategy.

The coupling approach between the Navier–Stokes solver and the Nernst–Planck–Poisson solver is based on a fixed point iteration strategy:

The projection of the pressure to the Poisson–Nernst–Planck solver just requires the evaluation of the pressure $p_{h}$ in the nodes of the triangulation.

According to (10), the projection of the velocity requires the integration of the reconstructed finite element solution $\Pi\vec{u}_{h}$ over interfaces between neighboring control volumes of the finite volume method. In the implementation, these integrals are calculated by quadrature rules. Sufficient accuracy of this step guarantees that the projected velocity is divergence free in the sense (11). For a detailed explanation of this algorithmically challenging step, see [27]. The projection operator can be assembled into a sparse matrix, allowing for computationally efficient repeated application of the projection operator during the fixed point iteration. As a consequence, in the case of electroneutral, inert transported species, the maximum principle for species concentrations is guaranteed, see [27] for more details. In combination with pressure robust finite element methods, this coupling approach was first applied to modeling of thin layer flow cells [31].

4 DC induced charge electroosmosis (ICEO) over an electrode with floating voltage

The discretization methods and the coupling strategy introduced above are implemented in the framework of the toolbox pdelib [32] that is developed at WIAS.

The solution of the Nernst–Planck–Poisson system is performed using Newton’s method. As the convergence of this method is guaranteed only in the vicinity of the solution, due to the strong nonlinearities, several measures need to be taken in order to obtain a solution. Time embedding is an approach to solve a stationary problem from a given initial value by solving a time dependent problem with increasing time steps by an implicit Euler method. Time step size control allows to guarantee that the difference between the solutions for two subsequent time steps is small enough to ensure the convergence of Newton’s method. Parameter ramping controls certain parameters (solvation number, surface charge) in the initial phase of the time evolution in such a way that the nonlinearities are easy to solve. For the flow part of the problem, the stationary Stokes problem is solved using a second order finite element method. Its velocity space consists of piecewise quadratic continuous vector fields enhanced with cell bubble functions and its pressure space consists of piecewise linear and discontinuous scalar fields [33]. For this method, an efficient divergence free reconstruction operator into the Raviart–Thomas finite element space of first order is constructed by standard interpolation [30, 34]. The method has been verified against classical Helmholtz–Smoluchowski theory and corresponding asymptotic expressions [35, 36, 37].

In order to demonstrate the capabilities of the model and the method, we simulate direct current (DC) induced charge electroosmosis (ICEO) over an electrode with a floating voltage. The effect of induced charge electroosmosis is based on a space charge region which appears at conducting surfaces immersed in an electric field. High electric conductivity keeps the potential of the conducting surface at a constant value. The external application of an electric field then must result in a potential gradient at the surface which triggers the formation of a space charge region. In this space charge region, electric forces act on the fluid, eventually setting the fluid in motion. For a more thorough discussion of this phenomenon including the alternating current case and other numerical approaches mainly based on the Helmholtz-Smoluchowski approximation for thin double layers see e.g. [38, 39, 40, 41].

We consider two large reservoirs that are separated by an impermeable wall with a narrow channel connecting the reservoirs. We assume that both reservoirs are filled with the same liquid electrolyte and that the ionic concentrations and the pressure are the same on both sides. For simplicity, we assume that the channel is of infinite width with planar parallel walls. Then, the computational domain $\Omega$ is restricted to the rectangular channel region in a cut plane, see Figure 2. On the boundary segment $\Gamma_{e}$ in the middle of the bottom wall there is a metal electrode which is not connected to any outer circuit.

The channel walls on the top and on the bottom side of $\Omega$ are impermeable, i.e. $\vec{N}_{i}\cdot\vec{n}=0$ and a no slip condition is imposed, i.e. $\vec{u}=0$. Except at the metal electrode $\Gamma_{e}$, the channel walls are electrically insulating and uncharged, i.e. $\varepsilon\nabla\psi\cdot\vec{n}=0$. The metal electrode $\Gamma_{e}$ is assumed to consist of an ideal metallic conductor with zero resistivity embedded in an insulating environment such that the electric potential within the metal is constant. The value of this constant is defined by the applied electric field. In this sense, the constant voltage at this electrode is floating. At the interfaces $\Gamma_{l}$ and $\Gamma_{r}$ to the reservoirs, a zero stress boundary condition $\nu\nabla\vec{u}\cdot\vec{n}=p\vec{n}$ for the flow is imposed which allows unimpeded motion of the fluid through the boundary of the simulation domain. A partial justification of this boundary condition consists in the fact that at these boundary we can ignore the electric forces due to local electroneutrality. Also, the same problem was simulated with a no-slip boundary condition $(\vec{u}=0)$ with similar results. Since the reservoir volume is considered large compared to the channel, we impose the rather crude approximating condition of constant prescribed species concentrations $c_{i}$ at the interfaces. The electric field is applied by imposing a potential difference between the electrolyte reservoirs. As an approximation, bias values of the same magnitude but opposite sign are applied at the reservoir boundaries $\Gamma_{l}$ and $\Gamma_{r}$. Due to symmetry, it then can be assumed that the floating potential value at the electrode on $\Gamma_{e}$ is zero: $\phi|_{\Gamma_{e}}=0$.

The potential difference between the reservoirs induces an electric current that is carried by (positively charged) cations moving from regions of higher potential to places with lower potential and (negatively charged) anions moving nearly in the opposite way. Without any charge accumulation, due to electroneutrality, there would be no influence on the velocity, and the flow velocity would be zero.

The electric field in lateral direction triggers the formation of a polarization boundary layer at the electrode. Due to symmetry, the charge density in this layer changes sign at the center of the electrode. As a consequence, there are electroosmotic forces of opposite sign acting on the fluid that lead to the creation of nanofluidic vortices. The stream plots in Figure 3 give a qualitative impression of the emerging flow and charge distribution for a 1M electrolytic solution. For the simulation, a series of discretization grids of increasing refinement level consisting of 465, 940,1875, 3690, 7596, 16383 nodes, respectively are used. The results of the modified and classical Nernst–Planck models can be compared. While one observes a significant difference in the width of the polarization boundary layer (due to the finite size limitations in the solvation based model), the flow pattern appears qualitatively very similar.

In order to obtain a more precise picture, Figures 4 and 5 exhibit profiles of the $y$-component of the flow. One observes good grid convergence behavior for both kinds of models, suggesting accurate numerical approximation of the respective solutions. The flow velocity for the solvation based model is considerably larger than for the classical Nernst–Planck model.

Figure 6 provides some insight into the influence of the solvation parameter $\kappa$ on the induced flow field. Increasing $\kappa$, a slight increase of the velocity is observed, once again consistent with the finite size effect.

Potential and concentration profiles in horizontal and in vertical direction are given in Figures 7 and 8. They demonstrate the strong differences between the classical Nernst-Planck model and the solvation based model with respect to the distribution of the electrostatic potential and the positive ions.

5 Conclusions

We presented a model for electro-osmotic flow at the nanoscale which includes finite ion size effects in a way that is consistent with first principles of non-equilibrium thermodynamics. We proposed a novel numerical solution approach which preserves important structural aspects of the continuous model. This method allows to simulate ion distributions in the vicinity of electrodes in such a way that ion concentrations in polarization boundary layers are not drastically overestimated and thus a meaningful coupling with the electro-osmotic flow due to the Coulomb force is possible. The capabilities of the numerical method have been demonstrated by a simulation of induced charge electroosmotic flow at a nanoscale electrode with floating potential. The simulation results for a 1M electrolyte show a considerable influence of finite ion size and solvation on electroosmotic flow velocity.

Nanoscale ion transport, electric field distribution and electrolyte flow are processes which need to be thoroughly understood in order to provide meaningful interpretations of experimental work in the emerging field of nanoscale electrochemistry. We hope to provide valuable contributions to future research efforts in this field by amending the model with more transported species and by incorporating electrochemical reactions at electrode surfaces into the model.

6 Acknowledgement

The research described in this paper has been supported by the German Federal Ministry of Education and Research Grant 03EK3027D (Network “Perspectives for Rechargeable Magnesium-Air batteries”) and the Einstein Foundation Berlin within the Matheon Project CH11 ”Sensing with Nanopores”.

7 Bibliography