Software Publications For Students Links About Author
PL
EN
For Students → PMiKNoM Seminars → Class 11 and 12
The content of this page is available under the Creative Commons Attribution-ShareAlike License. Valid XHTML 1.0 Transitional Valid CSS! Last modified: 26th of May 2020
PMiKNoM Seminars
Class 1
Class 2
Class 3
Class 4
Class 5
Class 6
Class 7 - Colloquium I
Class 8
Class 9
Class 10
Class 11 and 12
Class 13 - Colloquium II
Class 14 - The retake of the Colloquium
Class 15 - Assessment Colloquium
PMiKNoM Project
Sustainable Development
MSc and BSc theses
Print content

The Nernst-Planck-Poisson model

PMiKNoM seminars - classes 11 and 12

Introduction

Electro-diffusion is the motion of ions occuring as an effect of two processes:

Electro-diffusion processes are discussed in many fields of science, including semi-conductor science, analytical chemistry or neuro-physiology. Tool for the description of electro-diffusion is known as the Nernst-Planck-Poisson Model and can be applied to describe the processes occuring from nano to macro scale. For more details see:

NPP model - basic equations

The NPP model is an initial-boundary problem. Below the solution based on Dirichlet boundary conditions for concentrations and potential is presented.

Mass-balance equation

No reactions are considered, therefore the mass-balance equation takes the form:

$$ \frac{\partial c_i}{\partial t} = - \frac{\partial J_i}{\partial x} $$
where `c_i` represents the concentration of the `i`-th ion and `J_i` denotes its flux; `x` is space and `t` represents time.
Nernst-Planck equation

The Nernst-Planck equation for flux includes the description of diffusion and migration processes (as two components of the sum):

$$ J_i= - D_i \frac{\partial c_i}{\partial x} - \frac{F}{RT} z_i D_i c_i \frac{\partial \varphi}{\partial x} $$
where `D_i` represents the diffusion coefficient of the i-th ion and `z_i` denotes its charge; `\varphi` represents the electrical potential, `F` is the Faraday constant, `R` is the gas constant and `T` represents temperature.
Poisson equation

The Poisson equation relates the potential with the occuring deviation from electroneutrality:

$$ \frac{\partial^2 \varphi}{\partial x^2} = - \frac{F}{\varepsilon} \sum_i z_i c_i $$
where `\varepsilon` is the permittivity of the medium.
Initial and boundary conditions

Constant initial concentrations of all ions are assumed:

`c_i (x,0) = c_{i,ST} = const`

Dirichlet boundary conditions for potential are used, i.e. constant values of potential are assumed at the boundaries:

$$ \left\{ \begin{array}{l} \varphi(0,t) = \varphi_L = const \\ \varphi(l,t) = \varphi_R = const \\ \end{array} \right. $$
Dirichlet boundary conditions can be used also for concentrations:
$$ \left\{ \begin{array}{l} c_i(0,t) = c_{i,L} \cdot k_{in} / k_{out} = const \\ c_i(l,t) = c_{i,R} \cdot k_{in} / k_{out} = const \end{array} \right. $$

Alternatively Chang-Jaffe boundary conditions can be applied:

$$ \left\{ \begin{array}{l} J_i(0,t) = k_{in} \cdot c_{i,L} - k_{out} \cdot c_i(0,t)\\ J_i(l,t) = k_{out} \cdot c_i(l,t) - k_{in} \cdot c_{i,R} \end{array} \right. $$
where `k_{i n}` and `k_{out}` are the heterogeneous rate constants of `i`-th ion, entering and leaving the membrane, respectively.

NPP model - numerical solution

Numerical solution of the above equations is based on the Method of Lines (MoL) and on the explicit Euler method.

Method of Lines for concentrations

We divide the considered segment into `n+1` sub-segments, each having the length `h`, to obtain the discrete points `x_0` to `x_{n+1}`. Let's then define `c_{i,k} \equiv c_i(x_k)`, and define the additional points ` x_{k-1/2} ` and ` x_{k+1/2}` in the distance `h/2` from point `x_k`. In these points we will be establishing the fluxes (`J_{i,k-1/2}` and `J_{i,k+1/2}`).

Mass balance equation for the `i`-th component in `k`-th point of space is given as:

$$ \frac{dc_{i,k}}{dt} = - \frac{J_{i,k+1/2}-J_{i,k-1/2}}{h} $$

Nernst-Planck equation for fluxes `J_{i,k-1/2}` and `J_{i,k+1/2}` takes the form:

$$ \left\{ \begin{array}{l} J_{i,k-1/2}= - D_i \frac{d c_{i,k-1/2}}{dx} - \frac{F}{RT} z_i D_i c_{i,k-1/2} \frac{d \varphi_{k-1/2}}{dx} \\ J_{i,k+1/2}= - D_i \frac{d c_{i,k+1/2}}{dx} - \frac{F}{RT} z_i D_i c_{i,k+1/2} \frac{d \varphi_{k+1/2}}{dx} \end{array} \right. $$
Now the space derivatives, occuring in Nernst-Planck equation, can be substituted with finite differences:
$$ \left\{ \begin{array}{l} J_{i,k-1/2}= - D_i \frac{c_{i,k}-c_{i,k-1}}{h} - \frac{F}{RT} z_i D_i \frac{c_{i,k}+c_{i,k-1}}{2} \frac{\varphi_k -\varphi_{k-1}}{h} \\ J_{i,k+1/2}= - D_i \frac{c_{i,k+1}-c_{i,k}}{h} - \frac{F}{RT} z_i D_i \frac{c_{i,k+1}+c_{i,k}}{2} \frac{\varphi_{k+1} -\varphi_k}{h} \end{array} \right. $$

After substituting `J_{i,k-1/2}` and `J_{i,k+1/2}` into the mass balance equation, one obtains:

$$ \frac{dc_{i,k}}{dt} = \frac{D_i}{h^2} \left( c_{i,k+1}-2c_{i,k}+c_{i,k-1} \right) + \frac{F z_i D_i}{2RT h^2} \left( (c_{i,k+1}+c_{i,k})(\varphi_{k+1} -\varphi_k) - (c_{i,k}+c_{i,k-1})(\varphi_k -\varphi_{k-1}) \right). $$
Euler Method

If `c_{i,k}` denotes the concentration at time `t_j` and `c\text{_}n e w_{i,k}` denotes the the concentration at time `t_{j+1}` then:

$$ c\text{_}new_{i,k} = c_{i,k} + dt \cdot \frac{dc_{i,k}}{dt} $$

After inserting the time derivative, one obtains:

$$ c\text{_}new_{i,k} = c_{i,k} + \alpha_i \left( c_{i,k+1}-2c_{i,k}+c_{i,k-1} \right) + \beta_i \bigl( (c_{i,k+1}+c_{i,k})(\varphi_{k+1} -\varphi_k) - (c_{i,k}+c_{i,k-1})(\varphi_k -\varphi_{k-1}) \bigr) $$
where `\alpha_i = \frac{dt \cdot D_i}{h^2}` and `\beta_i = \frac{dt \cdot F z_i D_i}{2RT h^2}`
Method of Lines for Poisson equation

After substituting the space derivative with the finite difference, one obtains:

$$ \frac{\varphi_{k+1}-2\varphi_k+\varphi_{k-1}}{h^2} = - \frac{F}{\varepsilon} \sum_i z_i c_{i,k} $$

After then multiplying by `h^2`, the equation of the `i`-th element is obtained:

$$ \varphi_{k-1} - 2\varphi_k + \varphi_{k+1} = - \gamma \sum_i z_i c_{i,k} $$
where `\gamma = \frac{h^2 F}{\varepsilon}`

Knowing that ` \varphi_0 = \varphi_L = const` and `\varphi_{n+1} =\varphi_R =const`, we can substitute `k=1` and `k=n`, obtaining:

$$ -2\varphi_1 + \varphi_2 = - \varphi_L - \gamma \sum_i z_i c_{i,1} $$ $$ \varphi_{n-1} - 2\varphi_{n} = - \varphi_R - \gamma \sum_i z_i c_{i,n} $$

The set of these three equations represent the tri-diagonal matrix:

$$ \left[ \begin{array}{ccccccccc|c} -2 & 1 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & - \gamma \sum z_i c_{i,1} - \varphi_L \\ 1 & -2 & 1 & 0 & \cdots & 0 & 0 & 0 & 0 & - \gamma \sum z_i c_{i,2} \\ 0 & 1 & -2 & 1 & \cdots & 0 & 0 & 0 & 0 & - \gamma \sum z_i c_{i,3} \\ 0 & 0 & 1 & -2 & \cdots & 0 & 0 & 0 & 0 & - \gamma \sum z_i c_{i,4} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots\\ 0 & 0 & 0 & 0 & \cdots & -2 & 1 & 0 & 0 & - \gamma \sum z_i c_{i,n-3} \\ 0 & 0 & 0 & 0 & \cdots & 1 & -2 & 1 & 0 & - \gamma \sum z_i c_{i,n-2} \\ 0 & 0 & 0 & 0 & \cdots & 0 & 1 & -2 & 1 & - \gamma \sum z_i c_{i,n-1} \\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 1 & -2 & - \gamma \sum z_i c_{i,n} - \varphi_R \\ \end{array} \right] $$

NPP model - program's algorithm


Declarations - global constants and variables:
  1. Declare the global constants `l`, `t_{end}`, `p`, `n`, `m`, `F`, `R`, `T` and `\varepsilon`.
  2. Declare the two-dimensional arrays `c(1` to `p`, `0` to `n+1)` and `c__\text{new}(1` to `p`, `1` to `n)` as well as one-dimensional array `varphi(0` to `n+1)`.
  3. Declare the arrays `D(1` to `p)`, `z(1` to `p)`, `c_L(1` to `p)`, `c_{ST}(1` to `p)`, `c_R(1` to `p)`, and the variables `varphi_L` and `varphi_R`.
  4. Declare the variables `h`, `dt` and `gamma`, and arrays `alpha(1` to `p)` and `beta(1` to `p)`, and also variables `i`, `j`, `k`.

Procedure calculating the potential - Sub CalculatePotential():
  1. Declare the arrays `a(1` to `n-1)`, `b(1` to `n)`, `c c(1` to `n-1)`, `dd(1` to `n)` and variable `s`.
  2. Assign values on diagonal, i.e. for `k=1` to `n`   `dd(k)=-2`.
  3. Assign values above and below diagonal, i.e. for `k=1` to `n-1`:   `c c(k)=1`   and   `a(k)=1`.
  4. Assign values of free term b(k) using loop (for `k=1` to `n`):
    • reset variable `s`, i.e. `s = 0`,
    • calculate `\sum_i z_i c_{i,k}`, i.e. for `i=1` to `p`:   `s = s + z(i) \cdot c(i,k)`,
    • assign `b(k) = -\gamma \cdot s`.
  5. Modify `b` for the first and las row, i.e. `b(1)=b(1) - \varphi_L` and `b(n)=b(n) - \varphi_R`.
  6. Solve the tri-diagonal matrix (as it was done in class 4).
  7. Rewrite the result to the potential table, i.e. for `k=1` to `n`:   `\varphi (k)=b(k)`.

Main procedure - Sub NPP():
  1. Assign values to arrays `D`, `z`, `c_L`, `c_{ST}`, `c_R` and variables `varphi_L` and `varphi_R`.
  2. Calculate the step value `h=l/{n+1}`, the time-step value `dt=t_{end}/m` and the values of `alpha(i)`, `beta(i)` and `gamma`.
  3. Assign the initial values of concentrations `c(i,k)=c_{ST}(i)`. (embeded loop for `k=1` to `n` and for `i=1` to `p`).
  4. Assign the boundary values of concentrations, i.e. for `i=1` to `p`:   `c(i,0)=c_L(i)` and `c(i,n+1)=c_R(i)`.
  5. Assign the boundary values of potential: `varphi(0)=varphi_L` and `varphi(n+1)=varphi_R`.
  6. Call the procedure calculating potential, i.e. `Call` `Calcu latePotential`.
  7. Execute the embeded loop (for `j=1` to `m`), in which for every time-step:
    • program calculates the new values of concentration (for `k=1` to `n` and for `i=1` to `p`) and assign them to the array `c_\text{new}` (see equation 3.2.);
    • program switches the arrays, i.e. (for `k=1` to `n` and for `i=1` to `p`) `c(i,k)=c_\text{new}(i,k)`;
    • program calls the procedure calculating potential, i.e. `Call` `Calcu latePotential`.
  8. Save the results in the spreadsheet (in following columns: distance, concentrations of ions and potential).

NPP model - the obtained software