Electrokinetics part 1: advection-diffusion in 2D

Table of contents

1. Introduction
2. The advection-diffusion equation
3. Setting up the simulation
4. Comparison to the analytical solution
5. Numerical diffusion

This tutorial on the electrokinetics feature of ESPResSo is split into three parts. Each part is a self-contained notebook that treats one system:

• Part 1 (this notebook): advection-diffusion of a neutral tracer species in a prescribed two-dimensional flow field.
• Part 2 (electroosmotic flow): an electrolyte confined in a slit pore and driven by an external electric field.
• Part 3 (reactive flow): an irreversible bulk reaction in the turbulent wake behind an array of cylinders.

1. Introduction

In this tutorial we're looking at the electrokinetics feature of ESPResSo, which allows us to describe the motion of potentially charged chemical species solvated in a fluid on a continuum level. The governing equations for the solvent are known as the Poisson-Nernst-Planck equations, which is the combination of the electrostatic Poisson equation and the dynamics of the chemical species described by the Nernst-Planck equation. For the advection we solve the incompressible Navier-Stokes equation. The total set of equations is given by

$$ \begin{aligned} \partial_{t} n_{i} &= - \vec{\nabla} \cdot \vec{j}_{i} \\ \vec{j}_{i} &= - D_{i} \vec{\nabla} n_{i} - \frac{z_{i} e}{k_{B} T} n_{i} \vec{\nabla} \phi + n_{i} \vec{u} \\ \Delta \phi &= \frac{1}{4 \pi \varepsilon_{0} \varepsilon_{\mathrm{r}}} \sum_{i} z_{i} e n_{i} \\ \rho (\partial_{t} \vec{u} + (\vec{u} \cdot \vec{\nabla}) \vec{u}) &= - \vec{\nabla} p + \eta \Delta \vec{u} + \sum_{i} \frac{k_{B} T}{D_{i}} \vec{j}_{i} + \vec{f}_{\mathrm{ext}} \\ \vec{\nabla} \cdot \vec{u} &= 0, \end{aligned} $$

where $n_{i}$ denotes the ion density of species $i$, $\vec{j}_{i}$ the density flux, $D_{i}$ the diffusion coefficient, $z_{i}$ the valency, $e$ the elementary charge, $k_{B}$ the Boltzmann constant, $T$ the temperature, $\phi$ the electrostatic potential, $\varepsilon_{0}$ the vacuum permittivity, $\varepsilon_{\mathrm{r}}$ the relative permittivity, $\rho$ the fluid density, $\vec{u}$ the fluid velocity, $p$ the hydrostatic pressure, $\eta$ the dynamic viscosity, and $\vec{f}_{\mathrm{ext}}$ an external force density.

ESPResSo discretizes these equations on a lattice: the Navier-Stokes equation is solved with the lattice-Boltzmann (LB) method, while the Nernst-Planck equation for every species is solved with a finite-difference scheme on the same lattice. Both solvers are coupled through the advection term $n_{i} \vec{u}$ and through the friction force density $\sum_{i} k_{B} T \vec{j}_{i} / D_{i}$ that the species exert on the fluid.

2. The advection-diffusion equation

The system simulated in this part of the tutorial is the simple advection-diffusion of a drop of uncharged chemical species in a constant velocity field. To keep the computation time small, we restrict ourselves to a 2D problem, but the algorithm is also capable of solving the 3D advection-diffusion equation. Furthermore, we can also skip solving the electrostatic Poisson equation, since there are no charged species present. The equations we solve thus reduce to

$$ \partial_{t} n = D \Delta n - \vec{\nabla} \cdot (\vec{v} n). $$

The fundamental solution of this partial differential equation can be found analytically in the case of a constant velocity field $\vec{v}$ and a constant diffusion coefficient $D$. For a $d$-dimensional system, the solution of an initially infinitesimally small droplet at the origin can be written as

$$ n(\vec{x},t) = \frac{1}{(4 \pi D t)^{d/2}} \exp \left( - \frac{(\vec{x} - \vec{v} t)^2}{4 D t} \right). $$

The relative importance of the two transport mechanisms is measured by the Péclet number

$$ \mathrm{Pe} = \frac{\lVert \vec{v} \rVert L}{D}, $$

where $L$ is a characteristic length scale of the system. For $\mathrm{Pe} \gg 1$ transport is dominated by advection, for $\mathrm{Pe} \ll 1$ by diffusion. As we will see in the last section of this notebook, the advection-dominated regime is also the regime in which the discretization of the advection term becomes noticeable.

After importing the necessary packages, we start by defining the parameters of the simulation.

import espressomd
import espressomd.lb
import espressomd.electrokinetics

espressomd.assert_features(["WALBERLA"])

import numpy as np
import scipy.optimize

import matplotlib.pyplot as plt

plt.rcParams.update({'font.size': 14})

BOX_L = [80, 80, 1]
AGRID = 1.0
DIFFUSION_COEFFICIENT = 0.06
TAU = 1.0
EXT_FORCE_DENSITY = [0, 0, 0]
FLUID_DENSITY = 1.0
FLUID_VISCOSITY = 1.0
FLUID_VELOCITY = [0.04, 0.04, 0.0]
KT = 1.0

RUN_TIME = 400

3. Setting up the simulation

The third box dimension is a single lattice cell, which together with the periodic boundary conditions of ESPResSo turns the system into a quasi-two-dimensional one.

system = espressomd.System(box_l=BOX_L)
system.time_step = TAU
system.cell_system.skin = 0.4

print(f"Péclet number: {np.linalg.norm(FLUID_VELOCITY) * BOX_L[0] / DIFFUSION_COEFFICIENT:.0f}")

Péclet number: 75

The advecting flow field is provided by a lattice-Boltzmann fluid. Since we want a prescribed, constant velocity field, we initialize every lattice node with the same velocity and switch off the thermal fluctuations of the fluid.

Exercise 1

• Create an espressomd.lb.Lattice object called lattice with one ghost layer and lattice constant AGRID.
• Create an espressomd.lb.LBFluid called lbf on that lattice and register it with system.lb.
• Initialize the velocity of all lattice nodes to FLUID_VELOCITY.

Hints:

• Use the variables FLUID_DENSITY, FLUID_VISCOSITY, TAU and EXT_FORCE_DENSITY defined above.
• Set kT=0. to obtain a deterministic flow field. A thermalized fluid would superimpose random velocity fluctuations onto the constant advection velocity.
• Lattice nodes are addressed with a slice syntax, e.g. lbf[:, :, :].

# SOLUTION CELL
lattice = espressomd.lb.Lattice(n_ghost_layers=1, agrid=AGRID)
lbf = espressomd.lb.LBFluid(
    lattice=lattice, density=FLUID_DENSITY, kinematic_viscosity=FLUID_VISCOSITY,
    tau=TAU, ext_force_density=EXT_FORCE_DENSITY, kT=0.0, seed=42)
lbf[:, :, :].velocity = FLUID_VELOCITY
system.lb = lbf

To use the electrokinetics algorithm in ESPResSo, one needs to create an instance of the EKContainer object and pass it a time step tau and a Poisson solver solver. Since our species is uncharged, we don't need to solve the electrostatic Poisson equation, so we can use the placeholder class, which is called EKNone.

eksolver = espressomd.electrokinetics.EKNone(lattice=lattice, tau=TAU)
system.ekcontainer = espressomd.electrokinetics.EKContainer(tau=TAU, solver=eksolver)

Now we can add a diffusive species to the container and set its initial condition. To compare our simulation to the fundamental solution of the advection-diffusion equation, we need to approximate a delta-shaped droplet, which can be achieved by having a non-zero density only at the center of the domain. Individual lattice nodes of a species are accessed with the same slice syntax as the LB fluid, e.g. species[i, j, k].density.

Exercise 2

• Create an instance of espressomd.electrokinetics.EKSpecies called species and add it to the system with system.ekcontainer.add().
• Initialize the delta-shaped droplet by setting the density of the single lattice node at the center of the box to 1.0 / AGRID**3.

Hints:

• Use the variables DIFFUSION_COEFFICIENT, KT and TAU defined above.
• Enable both advection and friction_coupling.
• Make sure to construct the species with density=0, and disable electrostatics by setting valency to 0 as well.
• A single lattice node is addressed as species[i, j, k]; the center of the box in lattice units is [BOX_L[0] // 2, BOX_L[1] // 2, 0]. The species stores a density, so setting it to 1.0 / AGRID**3 places a total amount of solute of 1 on that node.

# SOLUTION CELL
species = espressomd.electrokinetics.EKSpecies(
    lattice=lattice, density=0.0, kT=KT,
    diffusion=DIFFUSION_COEFFICIENT, valency=0.0,
    advection=True, friction_coupling=True,
    ext_efield=[0., 0., 0.], tau=TAU)
system.ekcontainer.add(species)

# approximate a delta-shaped droplet carrying a total amount of solute of 1
species[BOX_L[0] // 2, BOX_L[1] // 2, 0].density = 1.0 / AGRID**3

Now everything is set and we can finally run the simulation by running the integrator. Afterwards we read the density field back off the species grid, again using slice access, and keep a copy for the analysis below.

system.integrator.run(RUN_TIME)

density_simulation = np.copy(species[:, :, 0].density)

4. Comparison to the analytical solution

For comparison, we evaluate the analytical solution on the lattice nodes. The node positions are measured relative to the center of the box, i.e. relative to the initial position of the droplet, and the analytical Gaussian is shifted by the distance $\vec{v} t$ travelled by the droplet.

def calc_gaussian(pos: np.ndarray, time: float, D: float):
    dim = pos.shape[-1]
    return (4 * np.pi * D * time)**(-dim / 2) * \
        np.exp(-np.sum(np.square(pos), axis=-1) / (4 * D * time))


# lattice node positions relative to the center of the box
xpos = (np.arange(BOX_L[0]) - BOX_L[0] // 2) * AGRID
ypos = (np.arange(BOX_L[1]) - BOX_L[1] // 2) * AGRID
positions_grid = np.stack(np.meshgrid(xpos, ypos, indexing="ij"), axis=-1)

# distance travelled by the droplet due to advection
drift = np.asarray(FLUID_VELOCITY[:2]) * RUN_TIME * TAU

analytic_density = calc_gaussian(pos=positions_grid - drift,
                                 time=RUN_TIME * TAU, D=DIFFUSION_COEFFICIENT)

fig, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(15, 7))

ax1.imshow(density_simulation, origin="lower", vmin=0, vmax=6e-3)
ax1.set_title("simulation")

imshow = ax2.imshow(analytic_density, origin="lower", vmin=0, vmax=6e-3)
ax2.set_title("analytic")
fig.colorbar(imshow, ax=[ax1, ax2], shrink=0.8)
plt.show()

Before looking at the shape of the droplet in detail, we check two properties that the advection-diffusion equation guarantees: the total amount of solute is conserved, and the center of mass of the distribution is transported with the flow velocity, $\vec{x}_{\mathrm{com}}(t) = \vec{v} t$. Both are simple reductions of the density field that we just read off the grid, so we compute them directly (the total amount is $\sum_{i} n_{i} a^{3}$, the center of mass $\sum_{i} n_{i} \vec{x}_{i} / \sum_{i} n_{i}$).

total_amount = np.sum(density_simulation) * AGRID**3
com = np.tensordot(density_simulation, positions_grid,
                   axes=([0, 1], [0, 1])) / np.sum(density_simulation)

print(f"total amount of solute: {total_amount:.6f} (initially 1.0)")
print(f"center of mass: simulation = {np.around(com, 3)}, analytic = {drift}")

total amount of solute: 1.000000 (initially 1.0)
center of mass: simulation = [15.849 15.849], analytic = [16. 16.]

The total amount of solute is conserved to machine precision, because the finite-difference scheme is written in terms of fluxes between neighboring lattice cells. The center of mass agrees with the analytical prediction to well below one lattice constant, i.e. the advection velocity is reproduced correctly.

To look at the shape of the droplet, we take a cut along the diagonal of the box. Because the flow velocity points along that diagonal, the peak of the droplet stays on it at all times.

values_diagonal = np.diagonal(density_simulation)
analytic_diagonal = np.diagonal(analytic_density)
positions_diagonal = np.arange(len(values_diagonal)) * np.sqrt(2) * AGRID


def gaussian_kernel(x, magnitude, mu, sigma):
    return magnitude * np.exp(-(x - mu)**2 / (2 * sigma**2))


popt, pcov = scipy.optimize.curve_fit(gaussian_kernel, positions_diagonal,
                                      analytic_diagonal, p0=[0.05, 70., 5.])
positions_analytic = np.concatenate([[positions_diagonal[0]],
                                     np.linspace(popt[1] - 5 * popt[2],
                                                 popt[1] + 5 * popt[2], 120),
                                     [positions_diagonal[-1]]])
values_analytic = gaussian_kernel(positions_analytic, *popt)

fig = plt.figure(figsize=(8, 5))
ax = fig.gca()
ax.plot(positions_diagonal, values_diagonal, "o", mfc="none", label="simulation")
ax.plot(positions_analytic, values_analytic, label="analytic")

ax.set_xlabel("position")
ax.set_ylabel("density")

plt.legend()
plt.show()

5. Numerical diffusion

From the plot one can see that the position of the density peak matches well. However, one also sees that the droplet in the simulation has spread more than it should. The reason is that the discretization used for the advection term introduces an artificial, additional diffusion to the system.

Let us quantify this numerical diffusion. Along the diagonal of the box the analytical solution reduces to a one-dimensional Gaussian: with $x = y = r / \sqrt{2}$ the exponent becomes $-(r - r_{0})^{2} / (4 D t)$, so the profile along the diagonal is a Gaussian of variance

$$ \sigma^{2} = 2 D t. $$

We therefore fit a Gaussian to the simulated diagonal profile and convert its width into an effective diffusion coefficient $D_{\mathrm{eff}} = \sigma^{2} / (2 t)$.

popt_sim, _ = scipy.optimize.curve_fit(gaussian_kernel, positions_diagonal,
                                       values_diagonal, p0=[0.05, 70., 5.])
diffusion_effective = popt_sim[2]**2 / (2. * RUN_TIME * TAU)

print(f"input diffusion coefficient:     {DIFFUSION_COEFFICIENT:.4f}")
print(f"effective diffusion coefficient: {diffusion_effective:.4f}")
print("numerical contribution:          "
      f"{100. * (diffusion_effective / DIFFUSION_COEFFICIENT - 1.):.1f} %")

input diffusion coefficient:     0.0600
effective diffusion coefficient: 0.0782
numerical contribution:          30.3 %

The effective diffusion coefficient is noticeably larger than the one we asked for. The excess is a property of the advection discretization and grows with the distance $\lVert \vec{v} \rVert \tau / a$ that the fluid travels per time step in units of the lattice constant. This is a fundamental limitation of the algorithm, which is why it cannot be applied to pure advection problems.

If you still have time, re-run this notebook with

• a smaller FLUID_VELOCITY, e.g. [0.01, 0.01, 0.0], and observe that the numerical contribution shrinks (keep in mind that the droplet then travels a shorter distance), or
• a smaller DIFFUSION_COEFFICIENT, i.e. a larger Péclet number, and observe that the relative error grows, because the numerical diffusion is independent of $D$.

In part 2 we will let the species carry a charge and let the flow field emerge self-consistently from the electrostatic forces instead of prescribing it.