# A Fast Stable Discretization of the

Constant–Convection–Diffusion–Reaction Equations

of Kinetic Capillary Electrophoresis (KCE)

###### Abstract

A discretization scheme is introduced for a set of convection–diffusion equations with a non-linear reaction term, where the convection velocity is constant for each reactant. This constancy allows a transformation to new spatial variables, which ensures the global stability of discretization. Convection–diffusion equations are notorious for their lack of stability, arising from the algebraic interaction of the convection and diffusion terms. Unexpectedly, our implemented numerical algorithm proves to be faster than computing exact solutions derived for a special case, while remaining reasonably accurate, as demonstrated in our runtime and error analysis.

MSC class: 65M08 (primary); 65M12, 35C05 (secondary).

Keywords: Convection–diffusion equations, multimesh, stable discretization.

mystylehidealllines=true,backgroundcolor=myyellow,innerleftmargin=0pt,innerrightmargin=0pt,innertopmargin=0pt,innerbottommargin=0pt,skipabove=0pt,skipbelow=0pt

###### Contents

## 1 Introduction

### 1.1 Aims and Overview

The aim of this paper is to develop a numerical algorithm for solving a general system of equations – the Constant–Convection–Diffusion–Reaction (CCDR) equations (Sec. 1.2) – applicable to the experimental models of both Kinetic Capillary Electrophoresis (KCE; Sec. 1.2) and KSEC (Sec. 4.2). Both are reversible binding reactions in a long and narrow tube: a capillary or a separation column. The CCDR equations model chemical reactions between any number of substances, with an arbitrary reaction term, convected in a fluid and propagated by a constant electric field, describing the spacetime evolution of reactant concentrations. Though the convection is constant in our model, this does not alone eliminate the instability of standard discretization schemes which convection–diffusion equations are notorious for, but it does allow a change of variables which induces stability (Sec. 2.1).

The presented mathematical endeavour is fundamentally important for two directions of our ongoing scientific research in kinetic separation: (1) the accurate computer simulation of experiments under appropriate conditions; (2) the resolution of an experimental–computational inverse problem for determining kinetic rate constants – both within the KCE framework. For these applications, a practical numerical solution method must be both stable and accurate, though for the inverse problem, a low runtime of this direct solver is also critical.

Previously, our simulations of KCE were generated with COMSOL, which employs a streamline upwind Petrov–Galerkin method to solve the general Nernst–Planck Equations (Secs. 1.2 and 2.1.1). This method is neither globally stable however, as evidenced by its output, nor fast enough for resolving the aforementioned inverse problem, described below. Furthermore, the use of a readily available black box simulation software, like COMSOL, has obvious limitations, such as the lack of flexibility for modification or extension to a broader framework, like an inverse problem. By their very generality, such software are likely to be suboptimal for a specific purpose. Thus the introduced numerical algorithm is not only an accurate simulation tool, but also an efficient direct solver serving a larger vision within the KCE framework, as follows.

The KCE system of convection–diffusion–reaction equations includes kinetic rate constants in the reaction term, which can be viewed as parameters that the solution functions of this system of partial differential equations depend on. The efficient and accurate generation of such parametric solutions to the direct problem is necessary for the evaluation of the error target function at each iteration of some optimization algorithm, which minimizes this target function to resolve the inverse problem of determining the parameters that induce a given target signal [30].

### 1.2 The Physical Model

The equations of our physical model are deduced from the Nernst–Planck Equations, which express the conservation of mass of ions in a fluid medium under the influence of an electric field along a single spatial direction, while accounting for convection, diffusion, and reaction between the ions.

Define (where ) as the spatiotemporal concentration of convected ions over points, as the velocities of their convection, as their diffusion, as a reaction term between the ions, as coefficients arising from various electromagnetic and thermodynamic constants [33], and as the electric fields influencing the motion of the ions, while denote with the Hadamard product of two vectors. To arrive at an intermediate form towards our equations, the Nernst–Planck equations [23, 33] reduce to the following vector PDE

with appropriate initial and boundary conditions that ensure the existence and uniqueness of a concentration vector solution.

For the experimental model of Kinetic Capillary Electrophoresis (KCE) [5, 21, 19, 15], the above equation reduces further due to certain features of the experimental setup. Specifically the convection velocities are constant , the diffusion coefficients are constant , and the electric field strength is constant as well. Thus a constant vector is introduced, and its components can be thought of as a case of constant convection in a convection–diffusion equation, suggesting the name Constant–Convection–Diffusion–Reaction (CCDR) equations. Therefore, they can be stated as

again with with appropriate initial and boundary conditions [25, 20].

In the formulation of the KCE model, there are three reactants and their concentration vector is denoted either as representing the experimental ligand, target, and complex, or as in the simplified case of MASKE [17]. The velocities and diffusion coefficients are denoted as and respectively. In the general KCE model, the reaction term is

Here are the rate constants of complex formation and dissociation respectively. In MASKE, the concentration of the second ion is assumed to be constant, so the second equation is omitted and the reaction term becomes . The point of this simplification is to enable us to find an exact solution [17]. Another simplification for the same reason, is to keep all three concentration components and simplify the reaction term only as , as in Sec. 2.2 [18].

### 1.3 Initial and Boundary Conditions

Since our goal is to give both a numerical and an exact solution to the above equations, the latter must be conveniently derivable from a Green function via function convolutions. So the initial and boundary conditions (IBC) of the CCDR equations are chosen to be idealizations of the initial concentration profiles measured in KCE experiments.

The initial conditions for the KCE equations specify the concentration profiles of the injected reactant plugs. A plug can be represented using various density functions . An identical substance amount is assumed within the plug – regardless of the density function used – to ensure that the amount of reacting molecules remains consistent for various injected concentration profiles. The total amount of molecules at concentration in a capillary of constant radius , at locations within cylinders of infinitesimal height , is given by the integral

So for two concentration functions to give the same amount of molecules, the areas under them must be equal. To be consistent with our earlier work [15], this amount is set at , which arises from the special case of rectangular concentration profiles of heights and plug length . To rescale a density function as another density, the proper transformation is , which retains a unit area under the curve. To ensure the standard amount of molecules, the area under the initial concentration profile must equal , implying the transformation .

This way our earlier table of initial and boundary conditions [15] can be generalized to arbitrary plug densities as follows.

Method | |||
---|---|---|---|

NECEEM | |||

cNECEEM | |||

SweepCE | |||

sSweepCE | |||

sSweepCEEM | |||

ECEEM | |||

ppKCE |

Note that denotes the location of the experimental detector at the end of the capillary, and denotes the Heaviside function

## 2 Solution of the Equations

### 2.1 Numerical Solution of the Complete Equations

#### 2.1.1 Discretization

As a reminder, the CCDR vector PDE equation is

The standard way to discretize the equations would be

yielding the vector iteration scheme

This common scheme is unstable due to the fact that the coefficients can be negative for a sparse mesh, causing oscillations to emerge in this recursion, commonly referred to as instability [31, 24, 12]. The main issue is that the Péclet number arising from the presence of both convection and diffusion (or the ratio ) is too large (convection-dominated case), requiring a very dense mesh for ensuring stability. This kind of instability can be partly resolved for convection–diffusion equations by weighting the diffusion term as in streamline upwind Petrov–Galerkin methods [13, 11, 1, 4]. Some programs have already been developed for modelling electrophoresis in other settings [28, 26, 8, 16, 27].

Instead of a general partial resolution, we may attempt to resolve the instability by virtue of constant convection in our equations. Petrov et al. and others [21, 9, 10] suggest improving stability by aligning the grid points with the directions of plug peak motion in spacetime, by choosing special individual stepsizes . Upon further thought, we can eliminate the convection term entirely by first changing variables for each line of the equations as

resulting in the transformed equations

Then the iterations will be executed separately for each concentration variable, and we interpolate between the meshes to compute the reaction term

We discretize the -planes as in order to align the grid points with the boundary conditions, which are now skewed in the transformed planes. So the iteration scheme becomes the following

The stability of the iteration is now ensured, since the terms in the earlier discretization are eliminated, so the condition for stability becomes which is easy to satisfy for a large Péclet number

So with , this implies a limit of at most grid points in time (and corresponding ). If , then there is no limitation for refinement.

#### 2.1.2 The Solver Algorithm

The algorithm presented below finds approximate values of the concentration vector solution of the CCDR equations with lines, over a given grid where , with an equidistant temporal discretization and arbitrary values in space.

A different -grid for each line of the equations is defined as where where (by choosing an equidistant spatial discretization ). The -grid for the transformed equations is defined as where (implying that for all ).

Denote the sought concentration values by , the values on the equidistant grid by , while the values on the transformed grids by . Note that , since by definition, for .

The reaction vector must be evaluated at the same spacetime grid points, so for each line of the equation we interpolate the concentration values on the other spatial grids onto the current -th grid , and denote the new values as (note that we can define ). The associated reaction values are denoted as

The algorithm contains nested for loops, first according to time, then by each reactant. Therefore, it begins by evaluating the initial conditions from some initial condition function . The left boundary condition is given by the function and the right boundary condition is of Neumann type implying that . Note that we must also extrapolate the transformed concentration for as

Lastly, denote the iteration coefficients (a convex combination) as

### 2.2 The Exact Solution for a Simplification of NECEEM

We plan to explicitly solve a simplified case of the CCDR vector PDE equation

with reaction function where , via the method of fundamental solutions as demonstrated earlier for the case of a rectangular plug [18]. Hereby the exact solutions are derived for a plug represented by a Gaussian initial condition, defined via the following density function

The NECEEM vector initial condition according to Sec. 1.3 is

while the boundary conditions remain unspecified for now, for the sake of simplicity.

Introducing some fundamental solutions will aid us in our derivation of an explicit solution. It can be shown by substitution that the equation

with the initial condition , is solved by

On the other hand, the fundamental solution satisfying

with a spatiotemporal Dirac delta function, is coincidentally the function .

Thus to find the third concentration component for some initial condition function , we must convolve it with the first fundamental solution above, resulting in

where is defined similarly as above, with parameters and .

In order to get the other two concentration components and , we break them up into an equilibrium term (the solution of the case above) and a dissipation term (the solution with right-hand side ), resulting in

Since the above formulas include several function convolutions (single and double integrals), it is desirable to simplify our derivation with the following well-known identity

Next, notice that for some and we have

Employing the above general solution formulas for this special initial condition, we get

With the above, we may now derive the other concentration components as follows

and are defined similarly.

Despite the extensive simplification of the exact solutions for the Gaussian case above – otherwise containing even more computationally expensive convolution integrals – we will show in Sec. 3.2 that somewhat surprisingly, the numerical method of the previous section is much faster according to our computational tests. It is not only faster, but is capable of solving the CCDR equations for any reaction function and initial-boundary conditions, with a reasonable error according to Sec. 3.3.

## 3 Computational Results

### 3.1 The Direct Solver Package

The numerical solver algorithm of Sec. 2.1.2 has been implemented in MATLAB for the KCE equations, though can be easily generalized to more reactants and any reaction mapping. The package is available on GitHub [29]. Currently the package contains only the direct solver for KCE, but the inverse solver built on it is also under development, serving our research described in the introduction.

The direct solver sub-package, contains an implementation of the numerical algorithm that is functional for all initial and boundary conditions, and plug types described in Sec. 1.3. Exact solutions have only been derived and implemented for the simplified NECEEM case in Sec. 2.2 with Gaussian and rectangular plug types [18], and for MASKE [17] with any plug type.

Therefore, in order to compare the computational runtimes and errors between the numerical and exact solutions, we are restricted to either of the above two simplified methods. We chose simplified NECEEM with a Gaussian plug for the analysis below, since its equations are closer to the full KCE equations than those of MASKE, consisting of only two lines. Furthermore, the solution formulas can be minimized computationally as much as possible, by eliminating expensive convolution integrals which arise for other density functions. This way we achieve a lower bound on the computational complexity of the exact case, which can be compared with that of the numerical one. Specifically, three single convolution integrals are eliminated for and two double convolution integrals are reduced to single ones for , as shown in Sec. 2.2.

The inverse problem of our research requires fast and accurate simulation of the solutions to the KCE equations for varying values, modeling experimental electropherogram signals [30]. Thus we analyze our direct solver in terms of both runtime and error below.

### 3.2 Runtime Analysis

In order to compare the runtime of our numerical solver (Sec. 2.1.2) to that of computing the exact solutions of simplified NECEEM (Sec. 2.2), we ran the solver for each temporal mesh size five times, and plotted the average of runtimes in Fig. 1. We are only interested in varying and not (the spatial mesh size), because the experimental signals are available only at the detector location , restricting the inverse problem to this single location in space. Test runs were done with an AMD A8-5550M processor, 8.00 GB of RAM, 64-bit Windows 10, and MATLAB R2016a, with the parameters

The runtimes seem to closely follow quadratic and linear relationships for the numerical and exact solvers, respectively, with high coefficients of determination . This is expected for the exact solver, which simply evaluates the formulas of Sec. 2.2. The quadratic polynomial for makes sense, since the algorithm of Sec. 2.1.2 must necessarily solve for all the nodes on the spacetime mesh to provide the values at the only relevant location . The transformed spatial meshes all depend on the temporal mesh linearly , implying that the spacetime mesh sizes must depend quadratically on .

The exact solution formulas for simplified NECEEM (Sec. 2.2) provide a standard of comparison for the numerical solution values, since the formulas for are exact, while the formulas for can be approximated accurately with some quadrature method (we used the integral function of MATLAB, with a RelTol value of ).

It is important to observe that the ratio of runtimes decreases for increasing , and follows a power law with an exponent of (though the ratio of a linear and a quadratic polynomial should asymptotically follow a hyperbola with an exponent of ). This implies that the gain in runtime using the numerical solver is much larger for lower mesh sizes, the difference ranging between 1-2 orders of magnitude below . For our end goal of an efficient inverse solver, a mesh size of is sufficient.

Thus we have demonstrated that surprisingly the numerical solver has a significantly lower runtime than even an optimally low runtime of the exact solution formulas. For this to be relevant, the error between the solutions must be reasonable.

### 3.3 Error Analysis

To show that the numerical solutions generated by Algorithm 1 are “reasonable”, we demonstrate that the -error from the exact values vanishes with an increasing temporal mesh size . Fig. 2 shows that the error decreases for increasing temporal mesh size according to similar power laws for the three concentration components. Taking the weighted average of the three exponents according to the coefficients of determination ( values) gives . Thus we conjecture that all three error components are inversely proportional to the temporal mesh size for large , i.e. precisely of order (Sec. 2.1.2).

The convergence of Algorithm 1 for increasing mesh size is self-evident from its design. Clearly, the error in the solution – arising from the derivative approximations, and the multimesh interpolation – vanishes as the spacetime mesh size increases. For an infinitesimally fine mesh, the error is necessarily zero. Thus the above comparison on Fig. 2 – between the exact and simulated solutions in the simplified case – is not intended to be a proof of convergence, but merely a rate analysis of the de facto vanishing of the error with respect to mesh size. It gives an idea of how “reasonable” this convergence is.

Nevertheless, this computational error analysis of the simplified system bears relevance to the original equations. The two differ in the reaction terms: (1) original: ; (2) simplified: , but the solutions differ negligibly. As illustrated on Fig. 3 for the simplified solution, the peaks are separated in a way that the following hold

Consequently, the following approximations hold

So we may conclude heuristically, that the above rate analysis of the error in the simplified case, is closely aligned with the original.

## 4 Concluding Remarks

### 4.1 Summary and Future Directions

A stable algorithm has been introduced for the efficient generation of accurate solutions to a set of convection–diffusion equations, with constant convection velocities. Stability has been demonstrated algebraically, due to the change of variables to new spatial meshes, for each line of this system of partial differential equations. Both efficiency and accuracy have been demonstrated as well.

The direct solver package corresponding to this work [29] may be extended in two potential directions: the implementation of new initial and boundary conditions, perhaps with a generic classification system; or the derivation of new exact solutions for either other simplifications of the original equations, or perhaps the explicit solution of the complete equations. According to our runtime analysis, however, the computation of closed form solutions is not likely to be more efficient, even in an extremely simple Gaussian case which gives a lower bound on computational complexity, implying the permanent significance of our solver.

The relevance of this algorithm to our further work on an experimental inverse problem has also been highlighted. Indeed, this is the direction our solver package is being extended in, now enabled by the fast computation of parametrized solutions with an error measured from an experimental signal, which is minimized via an optimization algorithm.

### 4.2 Applicability to Other Physical Models

The Kinetic Size-Exclusion Chromatography (KSEC) model [3, 2] can be written under certain experimental conditions in a form quite similar to the KCE equations [7, 6]. However, our program package in its current form [29] cannot handle the KSEC equations, mainly due to additional factors in the reaction term and the initial conditions, even in a simplified case [7]. Furthermore, due to their more special form, the KSEC equations can be reduced to an equivalent form, for which a more efficient solver should be developed.

### 4.3 Experimental Remarks

Both the KCE [15] and KSEC [3] equations describe kinetic separation methods, with concentration components as functions of time and space, where a reversible binding reaction is paralleled by the separation of reactants, convected at constant velocities in a capillary.

The presented discretization method was developed primarily for the KCE equations, which describe the migration of and interaction between three species: the ligand , the target , and the complex . The KCE equations do not incorporate a variety of other “background” species and reactions – such as components of the acid–base equilibration in the background electrolyte – which occur during electrophoresis of these three species.

In practice, the KCE method utilizes experimentally-measured velocities of the three species, and plug lengths corresponding to them. Such parameters are sufficient for the description of the relevant migration–reaction phenomena, only if the electrophoretic experiment is planned and carried out properly.

Most importantly, the concentrations of the analytes in KCE should be well-below the concentrations of the ions in the background electrolyte. In turn, in KSEC experiments, the equilibration between the pores and the free volume should be fast. Beneficially, these conditions are easily satisfied for both KCE and KSEC. Analyte concentrations are typically at least four orders of magnitude below the background electrolyte concentration, while the sub-nanometer size of pores allows fast equilibration between them and the free volume.

### 4.4 Acknowledgements

This work was supported by the Natural Sciences and Engineering Research Council of Canada (grant: CRDPJ 485321-15). We are grateful to Mirzo Kanoatov for noticing implementation bugs in the direct solver package, and for his suggestions towards improving its user-friendliness. We also appreciate the valuable suggestions made by the referees.

## References

- [1] R. E. Bank, J. F. Bürgler, W. Fichtner, and R. K. Smith. Some upwinding techniques for finite element approximations of convection-diffusion equations. Numerische Mathematik, 58(1):185–202, 1990.
- [2] J. Bao, S. M. Krylova, L. T. Cherney, J. Y. Le Blanc, P. Pribil, P. E. Johnson, D. J. Wilson, and S. N. Krylov. Pre-equilibration kinetic size-exclusion chromatography with mass spectrometry detection (peKSEC–MS) for label-free solution-based kinetic analysis of protein–small molecule interactions. Analyst, 140(4):990–994, 2015.
- [3] J. Bao, S. M. Krylova, L. T. Cherney, J. Y. LeBlanc, P. Pribil, P. E. Johnson, D. J. Wilson, and S. N. Krylov. Kinetic size-exclusion chromatography with mass spectrometry detection: An approach for solution-based label-free kinetic analysis of protein–small molecule interactions. Analytical chemistry, 86(20):10016–10020, 2014.
- [4] R. Becker and B. Vexler. Optimal control of the convection-diffusion equation using stabilized finite element methods. Numerische Mathematik, 106(3):349–367, 2007.
- [5] M. Berezovski and S. N. Krylov. Nonequilibrium capillary electrophoresis of equilibrium mixtures – a single experiment reveals equilibrium and kinetic parameters of protein-DNA interactions. Journal of the American Chemical Society, 124(46):13674–13675, 2002.
- [6] L. T. Cherney and S. N. Krylov. Slow-equilibration approximation in macroscopic approach to studying kinetics at equilibrium. Analytical chemistry, 83(4):1381–1387, 2011.
- [7] L. T. Cherney and S. N. Krylov. Slow-equilibration approximation in kinetic size exclusion chromatography. Analytical chemistry, 88(7):4063–4070, 2016.
- [8] P. Dutta, A. Beskok, and T. C. Warburton. Numerical simulation of mixed electroosmotic/pressure driven microflows. Numerical Heat Transfer: Part A: Applications, 41(2):131–148, 2002.
- [9] S. Ermakov, O. Mazhorova, and Y. Popov. Finite-difference algorithm for convection-diffusion equation applied to electrophoresis problem. Informatica, 3(2):173–197, 1992.
- [10] N. Fang and D. D. Chen. General approach to high-efficiency simulation of affinity capillary electrophoresis. Analytical chemistry, 77(3):840–847, 2005.
- [11] A. C. Galeão and E. G. D. Do Carmo. A consistent approximate upwind Petrov-Galerkin method for convection-dominated problems. Computer Methods in Applied Mechanics and Engineering, 68(1):83–95, 1988.
- [12] W. Hundsdorfer and J. G. Verwer. Numerical solution of time-dependent advection-diffusion-reaction equations, volume 33. Springer Science & Business Media, 2013.
- [13] C. Johnson. Numerical solution of partial differential equations by the finite element method. Courier Corporation, 2012.
- [14] M. Kanoatov, L. T. Cherney, and S. N. Krylov. Extracting kinetics from affinity capillary electrophoresis (ACE) data: A new blade for the old tool. Analytical chemistry, 86(2):1298–1305, 2014.
- [15] S. N. Krylov. Kinetic CE: Foundation for homogeneous kinetic affinity methods. Electrophoresis, 28(1-2):69–88, 2007.
- [16] L. Müllerová, P. Dubský, and B. Gaš. Twenty years of development of dual and multi-selector models in capillary electrophoresis: A review. Electrophoresis, 35(19):2688–2700, 2014.
- [17] V. Okhonin, M. V. Berezovski, and S. N. Krylov. MASKE: Macroscopic approach to studying kinetics at equilibrium. Journal of the American Chemical Society, 132(20):7062–7068, 2010.
- [18] V. Okhonin, S. M. Krylova, and S. N. Krylov. Nonequilibrium capillary electrophoresis of equilibrium mixtures, mathematical model. Analytical Chemistry, 76(5):1507–1512, 2004.
- [19] V. Okhonin, A. P. Petrov, M. Berezovski, and S. N. Krylov. Plug–plug kinetic capillary electrophoresis: Method for direct determination of rate constants of complex formation and dissociation. Analytical chemistry, 78(14):4803–4810, 2006.
- [20] O. Palusinski, A. Graham, R. Mosher, M. Bier, and D. Saville. Theory of electrophoretic separations. Part II: Construction of a numerical simulation scheme and its applications. AIChE journal, 32(2):215–223, 1986.
- [21] A. Petrov, V. Okhonin, M. Berezovski, and S. N. Krylov. Kinetic capillary electrophoresis (KCE): a conceptual platform for kinetic homogeneous affinity methods. Journal of the American Chemical Society, 127(48):17104–17110, 2005.
- [22] A. P. Petrov, L. T. Cherney, B. Dodgson, V. Okhonin, and S. N. Krylov. Separation-based approach to study dissociation kinetics of noncovalent DNA–multiple protein complexes. Journal of the American Chemical Society, 133(32):12486–12492, 2011.
- [23] R. F. Probstein. Physicochemical Hydrodynamics: An Introduction. John Wiley & Sons, New York, second edition, 1994.
- [24] H.-G. Roos, M. Stynes, and L. Tobiska. Robust numerical methods for singularly perturbed differential equations: convection-diffusion-reaction and flow problems, volume 24. Springer Science & Business Media, 2008.
- [25] D. Saville and O. Palusinski. Theory of electrophoretic separations. Part I: Formulation of a mathematical model. AIChE Journal, 32(2):207–214, 1986.
- [26] C. Schwer, B. Gaš, F. Lottspeich, and E. Kenndler. Computer simulation and experimental evaluation of on-column sample preconcentration in capillary zone electrophoresis by discontinuous buffer systems. Analytical Chemistry, 65(15):2108–2115, 1993.
- [27] S. Staal, M. Ungerer, A. Floris, H.-W. T. Brinke, R. Helmhout, M. Tellegen, K. Janssen, E. Karstens, C. van Arragon, S. Lenk, E. Staijen, J. Bartholomew, H. Krabbe, K. Movig, P. Dubský, A. van den Berg, and J. Eijkel. A versatile electrophoresis-based self-test platform. Electrophoresis, 36(5):712–721, 2015.
- [28] W. Thormann, M. C. Breadmore, J. Caslavska, and R. A. Mosher. Dynamic computer simulations of electrophoresis: a versatile research and teaching tool. Electrophoresis, 31(5):726–754, 2010.
- [29] J. Vass. KCE Solvers Package – Direct and Inverse Solver. GitHub/jzsfvss/KCESolvers, 2017.
- [30] J. Vass and S. N. Krylov. A computational resolution of the inverse problem of kinetic capillary electrophoresis (KCE). To be submitted, arXiv/1707.07852, 2017.
- [31] H. K. Versteeg and W. Malalasekera. An introduction to computational fluid dynamics: the finite volume method. Pearson Education, 2007.
- [32] Wikipedia. Finite difference method. Link, accessed 09-07-2017.
- [33] Wikipedia. Nernst-Planck equation. Link, accessed 09-06-2017.