Gouy-Chapman Debye layer

The Debye layer is the ionic concentration and potential distribution structure that appears on the surface of a charged electrode in contact with solvents in which are dissolved ionic species. Louis Georges Gouy and David Leonard Chapman at the beginning of the XX century proposed a model of the Debye layer resulting from the combined effect of its thermal diffusion and its electrostatic attraction or repulsion. In effect, in a stationary situation and assuming fluid at rest, the Poisson-Nernst-Planck equations are,


where ϕ is the electric potential and ci is the number of i-ions per volume. ωi and Zi are the i-ion mobility and valence. kB is the Boltzmann constant, e is the electron charge, ε the electrical permittivity and T the temperature.

The above equations, written in dimensionless form, reduces in the case of a fully dissolved binary system in a planar geometry to,


#include "grid/multigrid.h"
#include "diffusion.h"
#include "run.h"
#include "ehd/pnp.h"

#define Volt 1.0
#define DT 0.01

We assume a fully dissolved binary system labelling the positive ion as Cp and the counterion as Cm. The valence is one, (Z=1).

scalar φ[];
scalar Cp[], Cm[];
double dt;
int Z[2] = {1,-1};
scalar * sp = {Cp, Cm};

Ions are repelled by the electrode due to its positive volume conductivity while counterions are attracted (negative conductivity).

#if 1
const face vector kp[] = {1., 1.};
const face vector km[] = {-1., -1.};
vector * K = {kp, km};

On the left the charged planar electrode is set to a constant potential ϕ=1. The concentrations of the positive and negative ions depend exponentially on the voltage electrode.

φ[left] = dirichlet(Volt);
Cp[left]  = dirichlet (exp(-Volt));
Cm[left]  = dirichlet (exp(Volt));

In the bulk of the liquid, on the right boundary, the electrical potential is zero and the ion concentrations match the bulk concentration i.e

φ[right] = dirichlet (0.);
Cp[right]  = dirichlet (1.);
Cm[right]  = dirichlet (1.);

Initially, we set the ion concentration to their bulk values together with a linear decay of the electric potential ϕ.

event init (i = 0)
  foreach() {
    φ[] = Volt*(1.-x/5.);
    Cp[] = 1.0; 
    Cm[] = 1.0;
  boundary ({φ, Cp, Cm});

event integration (i++) {
  dt = dtnext (DT);

At each instant, the concentration of each species is updated taking into account the ohmic transport.

#if 1
  ohmic_flux (sp, Z, dt, K);
  ohmic_flux (sp, Z, dt); // fixme: this does not work yet

Then, the thermal diffusion is taken into account.

  for (scalar s in sp)
    diffusion (s, dt);

The electric potential ϕ has to be re-calculated since the net bulk charge has changed.

  scalar rhs[];
  foreach() {
    int i = 0;
    rhs[] = 0.;
    for (scalar s in sp)
      rhs[] -= Z[i++]*s[];
  poisson (φ, rhs);

event result (t = 3.5) {
    fprintf (stderr, "%g %g %g %g \n", x, φ[], Cp[], Cm[]);


We compare the numerical results (symbols) with the analytical solution (lines).

Profiles of electric potential and concentrations

Profiles of electric potential and concentrations

int main() {
  N = 32;
  L0 = 5;
  TOLERANCE = 1e-4;