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 socalled Poisson–Nernst–Planck (PNP) equations, which consist of the Poisson’s equation and the Nernst–Planck (NP) equations. Based on a meanfield 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 meanfield description, e.g., the steric effect, inhomogeneous dielectric effect, and ionion correlations [26, 23, 37, 55, 24, 28, 39, 30, 46, 51, 25, 52, 45, 33, 49].
In this work, we develop efficient and structurepreserving finite difference schemes for the PNP equations
(1) 
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 zeroflux boundary conditions possess several physically desired properties, including mass conservation, positivity preservation, and freeenergy dissipation, i.e.,
(2a)  
(2b)  
(2c) 
where the free energy , aside from some boundary contributions, is defined by
(3) 
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 (meanfield) 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 freeenergy 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 driftdiffusion 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 secondorder conservative, energy dissipative finite difference method was presented for the onedimensional 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:
(4) 
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 accuracypreserving 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 fixedpoint 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 NavierStokes 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 onedimensional 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 onedimensional 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 semiimplicit schemes that preserve positivity are much more computationally efficient, because larger timestepping sizes are allowed in temporal integration. In our recent work [11], we proposed efficient semiimplicit 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 semidiscrete form.
In this work, we develop implicit finite difference schemes for the multidimensional 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 centraldifferencing scheme based on harmonicmean 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 fixedpoint iteration method. Such a Netwon’s method can be employed to solve nonlinear systems resulting from other implicit discretization of the PNPtype 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 harmonicmean 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 structurepreserving 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, freeenergy 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 initialboundary value problem
(5) 
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](6) 
With initialboundary conditions given in (5), the property of freeenergy dissipation (2c) can be derived as well.
Theorem 2.1
The solution to the PNP equations (5) satisfies the energy dissipation law
(7) 

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 nonuniform 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
and
We denote by , , and the semidiscrete 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 dimensionbydimension manner.
3.2 Finite difference method in space
1) Spatial discretization of the Poisson’s equation
With given semidiscrete approximations , we discretize the Poisson’s equation with a central differencing stencil
(8) 
We employ central differencing stencils again to discretize the Dirichlet boundary conditions on by
(9) 
and the Neumann boundary conditions on by
(10) 
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
(11) 
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:
(12) 
3.3 Backward Euler method in time
To obtain numerical solutions of concentrations at different time steps, the semidiscrete 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 :
(17) 
which gives a fully implicit nonlinear system. Note that the zeroflux boundary conditions (16) have been used in (17). We further rewrite (17) in a matrix form as
(18) 
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 timedependent boundary input data. One efficient way to adjust the time step sizes has been proposed in [47, 29]:
(19) 
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.,
(20)  
(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
(22) We can verify that is an Mmatrix and in the elementwise 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
(23) 

From the proof of Theorem 4.2, we know that is an Mmatrix, and therefore exists and is a nonnegative 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 nonnegative, 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 freeenergy dissipation at full discrete level, when the boundary data is independent of time. The total discrete free energy (6) is approximated by
(24)  
Note that it is easy to verify that such an approximation is secondorder accurate in space.
Theorem 4.4
(Energy dissipation) The fully discrete free energy is nonincreasing for timeindependent boundary data in the sense that
(25) 
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
(26) where is between and , and is between and . Here the summation by parts have been used. By (26), we have
(27)
Comments
There are no comments yet.