Ion transport plays a fundamental role in many applications, such as electrochemical energy devices , electrokinetics in microfluidics , and transmembrane ion channels . 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., ). Long time behavior was studied in , 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 , 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 . A hybrid conservative scheme that uses adaptive grids was developed to solve the PNP equations on irregular domains . 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  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 , 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 . The resulting nonlinear discretization system was numerically solved by a fixed-point iteration method. Gao and He  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 .
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  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  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 , 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 
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 . 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
2) Spatial discretization of the NP equations
We first introduce a Slotboom reformulation of the Nernst–Planck equations:
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 :
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.
(Mass conservation)The fully implicit scheme (17) respect conservation of concentration, in the sense that the total concentration remains constant in time, i.e.,
(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 , 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 , and brief the corresponding proof with minor modifications.
The condition number of the coefficient matrix satisfies
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.
(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