Electrokinetics part 3: reaction in a turbulent flow

Table of contents

1. Introduction
2. Setting up the flow field
3. Setting up the reaction
4. Obstacles and sources
5. Running the simulation
6. Analysis

This is the third of the three parts of the electrokinetics tutorial:

• Part 1 (advection-diffusion): 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 (this notebook): an irreversible bulk reaction in the turbulent wake behind an array of cylinders.

This notebook is self-contained, but it assumes that you are familiar with the basic setup of an electrokinetics simulation as introduced in part 1.

1. Introduction

Besides advection and diffusion, the electrokinetics algorithm of ESPResSo can also treat chemical reactions between the species. The equation solved for every species $i$ is then an advection-diffusion-reaction equation,

$$ \partial_{t} n_{i} = - \vec{\nabla} \cdot \vec{j}_{i} + \nu_{i} \Gamma, \qquad \vec{j}_{i} = - D_{i} \vec{\nabla} n_{i} + n_{i} \vec{u}, $$

where $\nu_{i}$ is the stoichiometric coefficient of species $i$ and $\Gamma$ the reaction rate. Since none of the species in this part carries a charge, the electrostatic term drops out and the Poisson equation does not have to be solved.

To showcase the reaction feature, we simulate a simple reaction in a complex flow. For this, we choose a geometry of rigid cylinders. At large flow velocities, a Kármán vortex street, i.e., a repeating pattern of swirling vortices behind the obstacle, develops.

To this flow, we will add several species undergoing advection-diffusion, which is dominated by the downstream fluid flow in the channel. The reaction will be included as a bulk reaction, which means that the reaction can happen anywhere, the only requirement is that both species are present in the same lattice cell. When the reaction occurs, parts of the reactant species will turn into the product. How much of the species will transform within each time step is determined by the respective reaction rate and the overall structure of the reaction.

2. Setting up the flow field

We start by importing the necessary packages and defining the parameters.

import espressomd
import espressomd.lb
import espressomd.electrokinetics
import espressomd.shapes

espressomd.assert_features(["WALBERLA"])

import tqdm.auto as tqdm
import numpy as np

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.animation as animation

import tempfile
import base64

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

BOX_L = [80, 32, 1]
AGRID = 1.0
DIFFUSION_COEFFICIENT = 0.01
TAU = 0.03
EXT_FORCE_DENSITY = [0.6, 0, 0]
OBSTACLE_RADIUS = 6

DENSITY_FLUID = 0.5
VISCOSITY_KINEMATIC = 2.0
KT = 1.0

TOTAL_FRAMES = 100

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

The LB fluid is set up just as in part 1, so we give the code directly. There are two differences worth pointing out. First, the flow is driven along the $x$-direction by a constant body force density (ext_force_density). Second, the fluid is now thermalized: we pass a non-zero kT and register the fluid with the thermostat via system.thermostat.set_lb() (with gamma=0., since there are no particles). The thermal fluctuations break the symmetry of the flow around the obstacles and thereby trigger the vortex shedding.

lattice = espressomd.lb.Lattice(n_ghost_layers=1, agrid=AGRID)
lbf = espressomd.lb.LBFluid(
    lattice=lattice, density=DENSITY_FLUID, kinematic_viscosity=VISCOSITY_KINEMATIC,
    tau=TAU, ext_force_density=EXT_FORCE_DENSITY, kT=KT, seed=42)
system.lb = lbf
system.thermostat.set_lb(LB_fluid=lbf, seed=42, gamma=0.)

Since none of our species is charged, we can again use the placeholder Poisson solver EKNone.

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

3. Setting up the reaction

Now we can focus on the reactions. In this tutorial we choose the simple case of $A + B \rightarrow C$, which means that equal parts of the educt species $A$ and $B$ can turn into the product species $C$. ESPResSo distinguishes between educts and products by the sign of their respective stoichiometric coefficients, where educts have negative coefficients and products positive coefficients. Intuitively this can be understood that when a reaction happens, the density of the educts will decrease, hence the stoichiometric coefficient is negative.

The reaction rate constant $r$ is the rate at which the reaction happens. The order $O_i$ for a species $i$ specifies to which order the reaction depends on the density of that species. Positive orders mean that the reaction is faster the more density of this species is present, for negative orders the reaction slows down with higher density. In general, this process can be written as $\Gamma = r \left[ A \right]^{O_A} \left[ B \right]^{O_B} \left[ C \right]^{O_C}$, where $\Gamma$ is known as the reaction rate. This is sometimes also called the rate equation.

For our specific simulation this means that all stoichiometric coefficients are $-1$ for the educts and $+1$ for the product. We choose the order of the educts as $1$ and the order of the product as $0$. This means that the more amount of both educts is present, the more will react and the amount of product present won't have an influence.

REACTION_RATE_CONSTANT = 2.5
EDUCT_COEFFS = [-1, -1]
EDUCT_ORDERS = [1, 1]
PRODUCT_COEFFS = [1]
PRODUCT_ORDERS = [0]

We create each involved species and directly specify their boundary conditions for the domain boundaries. We set the initial density of the species to 0 and also add Dirichlet boundary conditions of zero density at both the inlet and the outlet of the system.

def create_species():
    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)
    # zero-density Dirichlet boundaries at the inlet and at the outlet
    species[0, :, :].density_boundary = espressomd.electrokinetics.DensityBoundary(0.0)
    species[-1, :, :].density_boundary = espressomd.electrokinetics.DensityBoundary(0.0)
    return species


educt_species = [create_species() for _ in EDUCT_COEFFS]
product_species = [create_species() for _ in PRODUCT_COEFFS]

Exercise 1

• Build a list reactants of EKReactant objects, one for every species that takes part in the reaction.
• Create an instance of EKBulkReaction using the reactants and activate the reaction by adding it to system.ekcontainer.reactions.

Hints:

• EKReactant takes the arguments ekspecies, stoech_coeff and order; use the lists EDUCT_COEFFS, EDUCT_ORDERS, PRODUCT_COEFFS and PRODUCT_ORDERS defined above.
• Educts and products are distinguished only by the sign of their stoichiometric coefficient, so a single flat list holds all of them.
• EKBulkReaction takes lattice, tau, the rate constant as coefficient and the list of reactants.

# SOLUTION CELL
reactants = []
for species, coeff, order in zip(educt_species, EDUCT_COEFFS, EDUCT_ORDERS):
    reactants.append(espressomd.electrokinetics.EKReactant(
        ekspecies=species, stoech_coeff=coeff, order=order))
for species, coeff, order in zip(product_species, PRODUCT_COEFFS, PRODUCT_ORDERS):
    reactants.append(espressomd.electrokinetics.EKReactant(
        ekspecies=species, stoech_coeff=coeff, order=order))

reaction = espressomd.electrokinetics.EKBulkReaction(
    reactants=reactants, coefficient=REACTION_RATE_CONSTANT,
    lattice=lattice, tau=TAU)

system.ekcontainer.reactions.add(reaction)

4. Obstacles and sources

The next thing to add to the system is the cylindrical obstacles, which act as the boundaries for the Kármán vortices to form. These are placed close to the inlet of the system and also act as impenetrable boundaries for the species. Since ESPResSo uses periodic boundary conditions, we need to add a total of three cylinders to the system, which will form two complete cylinders in the periodic system.

cylinder_centers = [
    [BOX_L[0] // 10, 0, 1],
    [BOX_L[0] // 10, BOX_L[1] // 2, 1],
    [BOX_L[0] // 10, BOX_L[1], 1],
]

shape_cylinder = []
for cylinder_center in cylinder_centers:
    shape_cylinder.append(espressomd.shapes.Cylinder(
        center=cylinder_center,
        axis=[0, 0, 1],
        length=BOX_L[2],
        radius=OBSTACLE_RADIUS,
    ))

We turn the cylinders into obstacles the same way the walls were set up in part 2, using add_boundary_from_shape on the fluid and on every species, so we give the code directly. Making a species impenetrable requires both a FluxBoundary with value [0, 0, 0] and a DensityBoundary with value 0..

for shape in shape_cylinder:
    lbf.add_boundary_from_shape(shape)
    for spec in (*educt_species, *product_species):
        spec.add_boundary_from_shape(
            shape, value=[0, 0, 0],
            boundary_type=espressomd.electrokinetics.FluxBoundary)
        spec.add_boundary_from_shape(
            shape, value=0.,
            boundary_type=espressomd.electrokinetics.DensityBoundary)

Up to this point there is no species present anywhere in the system and also no way for it to enter the system. Since the reaction is irreversible in our setup, we need to introduce some density of both the educt species to the system. For that we set two additional Dirichlet boundary conditions (sources) in the domain, where we fix the species' density to a constant, non-zero value. The sources are placed some distance apart along the transverse direction such that the reaction happens further downstream when the flow mixes the two species.

Exercise 2

• Place a source for the first educt at position source_x_pos and transverse position BOX_L[1] // 4, and a source for the second educt at the same downstream position but at 3 * (BOX_L[1] // 4). Both sources should be two lattice cells wide in the transverse direction.

Hints:

• A source is a DensityBoundary with a non-zero value; use the source_boundary object defined below.
• Assign it to the density_boundary property of the corresponding slice, e.g. species[i, j:k, :].density_boundary = source_boundary.

source_boundary = espressomd.electrokinetics.DensityBoundary(10.0)
source_x_pos = 1  # one lattice cell downstream of the inlet

# SOLUTION CELL
educt_species[0][source_x_pos, BOX_L[1] // 4 - 1:BOX_L[1] // 4 + 1, :].density_boundary = source_boundary
educt_species[1][source_x_pos, 3 * (BOX_L[1] // 4) - 1:3 * (BOX_L[1] // 4) + 1, :].density_boundary = source_boundary

5. Running the simulation

With this, the system is now finally complete and we can start the integration. To see the system evolve, we will render a movie from the timeseries of the system. For that we have to set up some helper functions for the plotting, which are beyond the scope of this tutorial.

VIDEO_TAG = """<video controls>
 <source src="data:video/x-m4v;base64,{0}" type="video/mp4">
 Your browser does not support the video tag.
</video>"""

# set ignore 'divide' and 'invalid' errors
# these occur when plotting the flowfield containing a zero velocity
np.seterr(divide='ignore', invalid='ignore')


def anim_to_html(anim):
    if not hasattr(anim, '_encoded_video'):
        with tempfile.NamedTemporaryFile(suffix='.mp4') as f:
            anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264'])
            with open(f.name, "rb") as g:
                video = g.read()
        anim._encoded_video = base64.b64encode(video).decode('ascii')
        plt.close(anim._fig)
    return VIDEO_TAG.format(anim._encoded_video)


animation.Animation._repr_html_ = anim_to_html

Nodes that belong to a boundary are excluded from the plots. The mask below is also used in the analysis section further down.

boundary_mask = lbf[:, :, 0].boundary != None

get_colormap = mpl.colormaps.get_cmap if hasattr(mpl.colormaps, "get_cmap") else mpl.cm.get_cmap
box_width = lattice.shape[1]
box_height = lattice.shape[0]

cmap = get_colormap("viridis").copy()
cmap.set_bad(color="gray")
cmap_quiver = get_colormap("binary").copy()
cmap_quiver.set_bad(color="gray")

# setup figure and prepare axes
fig = plt.figure(figsize=(9.8, 5.5))
imshow_kwargs = {"origin": "upper", "extent": (0, BOX_L[1], BOX_L[0], 0)}
gs = fig.add_gridspec(1, 4, wspace=0.1)
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1], sharey=ax1)
ax3 = plt.subplot(gs[2], sharey=ax1)
ax4 = plt.subplot(gs[3], sharey=ax1)
ax1.set_yticks(np.arange(0, BOX_L[0] + 1, 16))
for ax in (ax1, ax2, ax3, ax4):
    ax.set_xticks(np.arange(0, BOX_L[1] + 1, 16))

# set the background color for the quiver plot
bg_colors = np.copy(boundary_mask).astype(float)
bg_colors[boundary_mask] = np.nan
ax4.imshow(bg_colors, cmap=cmap_quiver, **imshow_kwargs)

for ax, title in zip(
    [ax1, ax2, ax3, ax4],
    ["educt 1", "educt 2", "product", "fluid velocity"]
):
    ax.set_title(title)
    ax.set_xlim((0, box_width))
    ax.set_ylim((0, box_height))

# create meshgrid for quiver plot subsampled by a factor 2
xs = np.arange(0, box_width, 2)
ys = np.arange(0, box_height, 2)
X, Y = np.meshgrid(xs, ys)

flowfield = lbf[:, :, 0].velocity[::2, ::2, :]
quiver = ax4.quiver(X + 1., Y + 1., flowfield[..., 1], flowfield[..., 0], scale=100)
fig.subplots_adjust(left=0.025, bottom=0.075, right=0.975, top=0.925, wspace=0.0125)

progress_bar = tqdm.tqdm(total=TOTAL_FRAMES)


def draw_frame(frame):
    system.integrator.run(50)

    flowfield = np.copy(lbf[:, :, 0].velocity)
    e1 = np.copy(educt_species[0][:, :, 0].density)
    e2 = np.copy(educt_species[1][:, :, 0].density)
    p = np.copy(product_species[0][:, :, 0].density)

    # apply the mask for the boundary
    e1[boundary_mask] = np.nan
    e2[boundary_mask] = np.nan
    p[boundary_mask] = np.nan
    flowfield[boundary_mask] = np.nan

    ax1.imshow(e1, cmap=cmap, vmin=0., vmax=source_boundary.density, **imshow_kwargs)
    ax2.imshow(e2, cmap=cmap, vmin=0., vmax=source_boundary.density, **imshow_kwargs)
    ax3.imshow(p, cmap=cmap, vmin=0., vmax=source_boundary.density, **imshow_kwargs)
    quiver.set_UVC((flowfield[::2, ::2, 1] + flowfield[1::2, 1::2, 1]) / 2.,
                   (flowfield[::2, ::2, 0] + flowfield[1::2, 1::2, 0]) / 2.)

    progress_bar.update()


animation.FuncAnimation(fig, draw_frame, frames=range(TOTAL_FRAMES), interval=300)

Looking at the movie of the species densities one can see that the fluid flow advects the educt species from their source locations past the cylinders into the system. Here, they start to mix and react, such that the product forms. The vortex street behind the obstacles enhances mixing through fluid turbulence. The density of the product then increases towards the outflow location of the channel, where it is deleted because of our zero-density boundary condition.

6. Analysis

To make the qualitative observations from the movie quantitative, we average the densities over the transverse direction and look at the resulting profiles along the flow direction. We also estimate the Reynolds number

$$ \mathrm{Re} = \frac{\bar{u}_{x} d}{\nu} $$

of the flow around the cylinders, where $\bar{u}_{x}$ is the mean downstream fluid velocity, $d = 2 R$ the cylinder diameter and $\nu$ the kinematic viscosity. For an isolated cylinder the wake becomes unsteady only above $\mathrm{Re} \approx 50$. Here the cylinders are spaced closely together — the gap between two of them is only a quarter of the transverse box length — and the resulting strong confinement destabilizes the wake already at a considerably lower Reynolds number.

The analysis itself is generic array post-processing of the density and velocity fields, so it is given here. The function mean_density_profile averages a species density over the transverse direction, excluding boundary nodes (boundary_mask is True on every boundary node; replacing those by np.nan and averaging with np.nanmean drops them). The mean downstream velocity follows the same pattern and gives the Reynolds number.

def mean_density_profile(species):
    density = np.copy(species[:, :, 0].density)
    return np.nanmean(np.where(boundary_mask, np.nan, density), axis=1)


positions = (np.arange(lattice.shape[0]) + 0.5) * AGRID

fig = plt.figure(figsize=(8, 5))
ax = fig.gca()
ax.plot(positions, mean_density_profile(educt_species[0]), label="educt 1")
ax.plot(positions, mean_density_profile(educt_species[1]), label="educt 2")
ax.plot(positions, mean_density_profile(product_species[0]), label="product")
ax.set_xlabel("downstream position")
ax.set_ylabel("transverse mean density")
ax.legend(loc="best")
plt.show()

velocity_mean = np.nanmean(
    np.where(boundary_mask, np.nan, lbf[:, :, 0].velocity[..., 0]))
reynolds_number = velocity_mean * 2. * OBSTACLE_RADIUS / VISCOSITY_KINEMATIC
print(f"mean downstream velocity: {velocity_mean:.4f}")
print(f"Reynolds number:          {reynolds_number:.1f}")

mean downstream velocity: 2.1712
Reynolds number:          13.0

The product profile rises from zero: right behind the sources the two educt streams are still separated, so no product can form, and the transverse mean of the product density then grows monotonically downstream as the vortices mix the two educts. The educt profiles, in contrast, are dominated by the strong density fluctuations of the vortex street — the reaction only consumes a modest fraction of the educts over the length of the channel. Note that these are instantaneous snapshots; averaging the profiles over many frames would give much smoother curves.

If you still have time, try the following variations:

• Increase REACTION_RATE_CONSTANT and observe how the product forms closer to the obstacles.
• Set EDUCT_ORDERS = [1, 0], i.e. make the reaction rate independent of the density of the second educt, and explain the resulting profiles.
• Reduce EXT_FORCE_DENSITY to lower the Reynolds number until the vortex street disappears, and observe how much less efficiently the two educt streams mix.