Comparison of drift-diffusion simulation with Gouy-Chapmann theory

In this note­book, di­lut­ed elec­trolyte be­tween metal­lic elec­trodes is sim­u­lat­ed.

The elec­tro­stat­ic po­ten­tial in elec­tri­cal dou­ble lay­er is plot­ted, and com­pared with an­a­lyt­i­cal Gouy-Chap­mann (GC) so­lu­tion giv­en by equa­tions:

\kappa=\sqrt{\frac{2c_0e^2}{\epsilon\epsilon_0k_BT}}

\psi=\frac{2k_BT}{e}\ *\ln \frac{e^{y_0/2}+1+(e^{y_0/2}-1)*\ e^{-\kappa x}}{e^{y_0/2}+1-(e^{y_0/2}-1)*e^{-\kappa x}}

y_0=\frac{e\psi_0}{k_BT}

with c_0 be­ing ion con­cen­tra­tion, \psi is the elec­tro­stat­ic po­ten­tial and oth­er sym­bols have usu­al mean­ing.

To en­sure ap­pli­ca­bil­i­ty of GC the­o­ry, low con­cen­tra­tion of ions is as­sumed. Dif­fer­ences be­tween the sim­u­la­tion and an­a­lyt­i­cal re­sults are caused by de­plet­ing elec­trolyte of ions.

%matplotlib inline
import matplotlib.pylab as plt
from oedes import bdf1adapt, models, testing, progressbar, init_notebook
init_notebook()
from oedes.fvm import mesh1d
import numpy as np
import scipy.constants
def run(mesh, c, appl_v):
    model = models.BaseModel()
    models.add_ions(model, mesh, zc=1,za=-1)
    model.setUp()
    params = {'T': 300.,
              'electrode0.voltage': 0,
              'electrode1.voltage': 0,
              'electrode0.workfunction': 0,
              'electrode1.workfunction': 0,
              'cation.mu': 1e-11,
              'anion.mu': 1e-11,
              'epsilon_r': 3.}
    xinit = model.X.copy()
    for eq in model.species:
        xinit[eq.idx] = c
    params['electrode0.voltage'] = appl_v
    for t, x, xt, out in bdf1adapt(model, xinit, params, 0., t1=1e10, dt=1e-9):
        o = out()
    return mesh, o
def edl_solution(appl_v, x, c, T=300., epsi_r=3.):
    kT_q = scipy.constants.Boltzmann * T / scipy.constants.elementary_charge
    K = np.sqrt(2 * c * scipy.constants.elementary_charge /
                (kT_q * scipy.constants.epsilon_0 * epsi_r))
    y0 = appl_v * 0.5 / kT_q
    t = np.exp(y0 / 2.)
    u = (t - 1.) * np.exp(-K * x)
    v = t + 1
    psi = 2 * kT_q * np.log((v + u) / (v - u))
    return psi
def test(appl_v):
    c = 1e24
    mesh, o = run(mesh1d(1e-6, dx_boundary=1e-14), c, appl_v)
    x = mesh.cells['center']
    psi = edl_solution(appl_v, x, c)
    testing.store(o['potential'] - appl_v / 2., rtol=1e-7, atol=1e-7)
    plt.plot(mesh.cells['center'], o['potential'] - appl_v / 2., '-')
    plt.plot(x, psi, '--')
    plt.xscale('log')
for v in progressbar([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7],'V'):
    test(v)
plt.ylim([0, 0.4])
plt.xlim([1e-14, 0.5e-6])
plt.xlabel('x')
plt.ylabel('potential')