Structure-Preserving and Efficient Numerical Methods for Ion Transport

by   Jie Ding, et al.

Ion transport, often described by the Poisson–Nernst–Planck (PNP) equations, is ubiquitous in electrochemical devices and many biological processes of significance. In this work, we develop conservative, positivity-preserving, energy dissipating, and implicit finite difference schemes for solving the multi-dimensional PNP equations with multiple ionic species. A central-differencing discretization based on harmonic-mean approximations is employed for the Nernst–Planck (NP) equations. The backward Euler discretization in time is employed to derive a fully implicit nonlinear system, which is efficiently solved by a newly proposed Newton's method. The improved computational efficiency of the Newton's method originates from the usage of the electrostatic potential as the iteration variable, rather than the unknowns of the nonlinear system that involves both the potential and concentration of multiple ionic species. Numerical analysis proves that the numerical schemes respect three desired analytical properties (conservation, positivity preserving, and energy dissipation) fully discretely. Based on advantages brought by the harmonic-mean approximations, we are able to establish estimate on the upper bound of condition numbers of coefficient matrices in linear systems that are solved iteratively. The solvability and stability of the linearized problem in the Newton's method are rigorously established as well. Numerical tests are performed to confirm the anticipated numerical accuracy, computational efficiency, and structure-preserving properties of the developed schemes. Adaptive time stepping is implemented for further efficiency improvement. Finally, the proposed numerical approaches are applied to characterize ion transport subject to a sinusoidal applied potential.



There are no comments yet.


page 21


Convergence Analysis of Structure-Preserving Numerical Methods Based on Slotboom Transformation for the Poisson–Nernst–Planck Equations

The analysis of structure-preserving numerical methods for the Poisson–N...

On higher order passivity preserving schemes for nonlinear Maxwell's equations

We present two strategies for designing passivity preserving higher orde...

An Assessment of Solvers for Algebraically Stabilized Discretizations of Convection-Diffusion-Reaction Equations

We consider flux-corrected finite element discretizations of 3D convecti...

Two novel classes of energy-preserving numerical approximations for the sine-Gordon equation with Neumann boundary conditions

We develop two novel classes of energy-preserving algorithms for the sin...

Preconditioning Richards Equations: spectral analysis and parallel solution at very large scale

We consider here a cell-centered finite difference approximation of the ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Ion transport plays a fundamental role in many applications, such as electrochemical energy devices [4], electrokinetics in microfluidics [48], and transmembrane ion channels [53]. It is often described by the so-called Poisson–Nernst–Planck (PNP) equations, which consist of the Poisson’s equation and the Nernst–Planck (NP) equations. Based on a mean-field approximation, the NP equations describe the diffusion of ions in the gradient of the electrostatic potential. The Poisson’s equation determines the electrostatic potential with the charge density arising from diffusing ions. Recently, there has been growing interests in incorporating effects that are beyond the mean-field description, e.g., the steric effect, inhomogeneous dielectric effect, and ion-ion correlations [26, 23, 37, 55, 24, 28, 39, 30, 46, 51, 25, 52, 45, 33, 49].

In this work, we develop efficient and structure-preserving finite difference schemes for the PNP equations


Here is the ion concentration for the -th species, is the valence of the -th ionic species, is a coefficient arising from nondimensionalization, is the electrostatic potential, and is the fixed charge density.

The analytical solutions to (1) with zero-flux boundary conditions possess several physically desired properties, including mass conservation, positivity preservation, and free-energy dissipation, i.e.,


where the free energy , aside from some boundary contributions, is defined by


The free energy contains both an entropic contribution and an electrostatic energy: is the entropy related to the Brownian motion of each ion species, and is the (mean-field) electrostatic energy of the Coulomb interaction between charged ions. The concentrations are expected to converge to an equilibrium solution in a closed system regardless of how initial data are distributed.

These nice mathematical features are crucial for the analytical study of the PNP equations. For instance, by an energy estimate with the control of the free-energy dissipation, the solution is shown to converge to the thermal equilibrium state as time becomes large, if the boundary conditions are in thermal equilibrium (see, e.g., [16]). Long time behavior was studied in [6], and further in [2, 5] with refined convergence rates. Results for the drift-diffusion model, i.e., the PNP equations in the semiconductor literature, with regarding global existence, uniqueness, and asymptotic behavior in the case of different boundary conditions have been established in the works [42, 17, 13, 12].

The PNP equations can hardly be solved analytically due to the nonlinear coupling of the electrostatic potential and ionic concentrations. Much effort has been devoted to the development of numerical methods in various applications [7, 43, 36, 54, 15, 8, 31, 41, 35, 50, 40, 14, 18, 32, 9, 49, 44, 10]. The existing algorithms range from finite difference to finite elements in both one dimension and high dimensions. Among these developed schemes, several attempts are made to design desirable numerical schemes that respect the nice properties (2) of analytical solutions. In [15], for instance, a second-order conservative, energy dissipative finite difference method was presented for the one-dimensional PNP equations. A delicate temporal discretization scheme was designed to preserve energy dynamics in [14]. A hybrid conservative scheme that uses adaptive grids was developed to solve the PNP equations on irregular domains [41]. A type of finite difference schemes has been developed using the Slotboom transformation of the Nernst–Planck (NP) equations:


where are the Slotboom variables [38, 36, 39, 32, 33]. By using the Slotboom variables, Liu and Wang [31] developed a free energy satisfying finite difference scheme that preserves those three properties with rigorous proof in one dimension. A free energy satisfying discontinuous Galerkin method was also developed by the same authors in [32], in which the positivity of numerical solutions was not proved but enforced by an accuracy-preserving limiter. An implicit finite difference scheme was developed to solve the PNP equations with properties of positivity preservation and energy dissipation [22]. The resulting nonlinear discretization system was numerically solved by a fixed-point iteration method. Gao and He [18] proposed a linearized convergent finite element scheme that conserves total concentration and preserves the electric energy. Rigorous error analysis of finite element type methods for the PNP equations has been studied in [50, 19]. A finite element discretization that can enforce positivity of numerical solutions was developed for the PNP equations, as well as the PNP equations coupling with the incompressible Navier-Stokes equations [40].

Although some progress has been made on the development of numerical methods that ensure the desired properties, it is still desirable to have computationally efficient and robust finite difference schemes that incorporate all three desired properties together, especially in high dimensions. Among three properties (2), the preservation of positivity is crucial to the validity of a numerical solution and is in particular hard to achieve. For instance, the positivity has been proved in [15] for the one-dimensional PNP equations, under assumptions that the gradient of the electrostatic potential is bounded and mesh step sizes satisfy certain constraint conditions. The numerical scheme developed in [31] has been proved to preserve positivity of the numerical solution but in one-dimensional case. Furthermore, one constraint on a mesh ratio needs to be satisfied to ensure positivity, due to the explicit nature of the scheme. From a practical point of view, implicit or semi-implicit schemes that preserve positivity are much more computationally efficient, because larger time-stepping sizes are allowed in temporal integration. In our recent work [11], we proposed efficient semi-implicit schemes that respect the desired properties. While the scheme allows relatively large time steps and is successful in positivity persevering and mass conservation, the energy dissipation is only proved in a semi-discrete form.

In this work, we develop implicit finite difference schemes for the multi-dimensional PNP equations with multiple ionic species that respect conservation, positivity preserving, and energy dissipation at fully discrete level. The NP equations reformulated in the Slotboom variables are spatially discretized by a central-differencing scheme based on harmonic-mean approximations. The backward Euler method in time is employed to derive a nonlinear coupled system, which is efficiently solved by a newly proposed Newton’s method. The improved efficiency of the Newton’s method is achieved via using the electrostatic potential as the iteration variable, rather than the unknowns of the nonlinear system that involves both the potential and concentrations of multiple ionic species. The advantages of the proposed Newton’s method in saving memory and computational efficiency become more significant when the number of ionic species gets larger. In addition, numerical simulations demonstrate that the Newton’s method requires appreciably fewer iteration steps and less computational time, in comparison with a typical fixed-point iteration method. Such a Netwon’s method can be employed to solve nonlinear systems resulting from other implicit discretization of the PNP-type equations.

We perform detailed numerical analysis to prove that the numerical schemes respect three desired analytical properties fully discretely. Thanks to the advantages brought by the harmonic-mean approximations, we are able to establish upper bounds on condition numbers of coefficient matrices of linear systems resulting from both the discretization of the NP equations, and solvability and stability of the linearized problem in the Newton’s method. The linear systems are efficiently solved using iterative methods with preconditioners. Numerical simulations are presented to demonstrate that the developed schemes have expected numerical accuracy, computational efficiency, and structure-preserving properties. An adaptive time stepping strategy is employed to achieve further improvement in computational efficiency. The benefit of the adaptive time stepping is demonstrated in the application of the proposed numerical approaches to probing ion transport in response to an alternating applied potential. Finally, the developed numerical approaches are applied to understand charge dynamics in electrolytes between two parallel electrodes exposed to sinusoidal applied potentials. The impact of frequency of the sinusoidal applied potentials is extensively investigated in numerical simulations.

The rest of this paper is organized as follows. In Section §3, we start with details of our settings, and present our implicit finite difference method in both spatial and temporal discretization. In Section §4 we prove the main properties at fully discrete level, including conservation, free-energy dissipation, positivity preservation. Section §5 is devoted to numerical examples, including accuracy test, charge dynamics, and adaptive time stepping. Finally, we conclude in Section §6.

2 The PNP equations

We consider the initial-boundary value problem


Here is a bounded domain, n

is a unit exterior normal vector on the boundary

, are initial concentration distributions. To be general, we here consider both Dirichlet and Neumann boundary conditions for the electrostatic potential, i.e., is a given electrostatic potential defined on the Dirichlet boundary , and is the surface charge density defined on the Neumann boundary with and . The corresponding total free energy with boundary contributions is given by [34]


With initial-boundary conditions given in (5), the property of free-energy dissipation (2c) can be derived as well.

Theorem 2.1

The solution to the PNP equations (5) satisfies the energy dissipation law

  • Taking a derivative with respect to time, we have by integration by parts that

    It follows from the Poisson’s equation that

    Substituting into the above equation completes the proof.       

One can observe from (7) that the total free energy of the system is dissipating if the boundary data are independent of time.

3 Numerical Method

3.1 Discretization and Notations

The spatial discretization is similar to our previous work [11]. For completeness of the presentation, we briefly introduce notations and recall the discretization scheme. For simplicity, we consider a 2D rectangular computational domain with boundaries

The computational domain is covered by the non-uniform grid points with

where and are the number of grid points along each dimension. We also introduce grid points with integer indices:

The grid spacings are given by


We denote by , , and the semi-discrete approximations of , , and , respectively. Define discrete operators

Discrete operators and can be analogously defined. Also, we introduce

We now present our implicit finite difference method by discretiztion in space and time separately. Note that, although our discretization is presented for a 2D case, it can be readily extended to three dimensions in a dimension-by-dimension manner.

3.2 Finite difference method in space

1) Spatial discretization of the Poisson’s equation

With given semi-discrete approximations , we discretize the Poisson’s equation with a central differencing stencil


We employ central differencing stencils again to discretize the Dirichlet boundary conditions on by


and the Neumann boundary conditions on by


Notice that the definition of the boundary data and have been extended to the whole computational domain. In numerical implementation, the ghost points outside are eliminated by coupling the discretization scheme (8) and boundary discretization (9) and (10). The coupled difference equations can be written in a matrix form


Here is the coefficient matrix, the column vector results from the boundary conditions (9) and (10), and , , and are vectors with components being , , and , respectively.

2) Spatial discretization of the NP equations

We first introduce a Slotboom reformulation of the Nernst–Planck equations:


where and , which are the Slotboom variables [36, 32].

It follows from central-differencing discretization of Eq. (12) at that


where the flux is approximated by


Harmonic-mean approximations are proposed in [11] to approximate the exponential terms at half-grid points:


The zero-flux boundary conditions are discretized by


3.3 Backward Euler method in time

To obtain numerical solutions of concentrations at different time steps, the semi-discrete scheme can be integrated with various ODE solvers implicitly. With a nonuniform time step size and , , , and are used to denote the numerical approximations of , , and , respectively.

We employ the backward Euler discretization in :


which gives a fully implicit nonlinear system. Note that the zero-flux boundary conditions (16) have been used in (17). We further rewrite (17) in a matrix form as


Here is a square matrix dependent on and , is a diagonal matrix given by

and and are column vectors with components being and , respectively.

An adaptive time stepping strategy is often useful in speeding up simulations of charged systems that have time-dependent boundary input data. One efficient way to adjust the time step sizes has been proposed in [47, 29]:


where is a given constant, is the free energy defined by (6). In our numerical implementation, the temporal derivative of the free energy is approximated by a difference quotient. The time steps and dictate the lower and upper bounds of the adaptive time steps, respectively, i.e., .

4 Properties of numerical solutions

The numerical solution has several important properties, as stated in the following theorems.

Theorem 4.1

(Mass conservation)The fully implicit scheme (17) respect conservation of concentration, in the sense that the total concentration remains constant in time, i.e.,

  • It follows from (13) that

    where we have used the zero-flux boundary conditions (16) in the last step. Similarly, summing both sides of (17) over gives

    Incorporating the zero-flux boundary conditions (16) at time for leads to (21).       

Theorem 4.2

(Positivity preserving)The numerical solutions computed from the backward Euler scheme (17) remain positive in time, i.e., if , then

  • Analogous to the derivation in our previous work [11], we can show that


    We can verify that is an M-matrix and in the element-wise sense. Thus, if .      

The concentrations are obtained by solving the linear system (18) iteratively. For iterative methods, it is desirable to establish estimates on the condition number of the coefficient matrix. We now recall the estimate on the condition number of the coefficient matrix in the work [11], and brief the corresponding proof with minor modifications.

Theorem 4.3

The condition number of the coefficient matrix satisfies

  • From the proof of Theorem 4.2, we know that is an M-matrix, and therefore exists and is a non-negative matrix. As shown in [11], we can further prove that

    where . This implies that each column sum of is less or equal to . Since each element in is non-negative, we obtain . Also, it follows from (22) that

    This completes the proof by the definition of the -norm condition number.       

Our numerical scheme also respects the property of free-energy dissipation at full discrete level, when the boundary data is independent of time. The total discrete free energy (6) is approximated by


Note that it is easy to verify that such an approximation is second-order accurate in space.

Theorem 4.4

(Energy dissipation) The fully discrete free energy is non-increasing for time-independent boundary data in the sense that


where is a number between and , and is a number between and .

  • The fully discrete NP equations are given by

    Multiplying both sides by and summing over indices lead to


    where is between and , and is between and . Here the summation by parts have been used. By (26), we have