J. Comput.
Phys.
Journal of Computational Physics 00 (2017) 1–33

arXiv:1708.08741v1 [] 29 Aug 2017

Coupled Multiphysics Simulations of Charged Particle
Electrophoresis for Massively Parallel Supercomputers
Dominik Bartuschata,∗, Ulrich Rüdea,b
a Lehrstuhl

für Systemsimulation, Friedrich-Alexander Universität Erlangen-Nürnberg, Cauerstrasse 11, 91058 Erlangen,
Germany
b Parallel Algorithms Group, CERFACS, 42 Avenue Gaspard Coriolis, 31057 Toulouse, France

Abstract
The article deals with the multiphysics simulation of electrokinetic flows. When charged particles are immersed in a
fluid and are additionally subjected to electric fields, this results in a complex coupling of several physical phenomena.
In a direct numerical simulation, the dynamics of moving and geometrically resolved particles, the hydrodynamics of
the fluid, and the electric field must be suitably resolved and their coupling must be realized algorithmically. Here
the two-relaxation-time variant of the lattice Boltzmann method is employed together with a momentum-exchange
coupling to the particulate phase. For the electric field that varies in time according to the particle trajectories,
a quasistatic continuum model and its discretization with finite volumes is chosen. This field is coupled to the
particulate phase in the form of an acceleration due to electrostatic forces and conversely via the respective charges
as boundary conditions for the electric potential equation. The electric field is also coupled to the fluid phase by
modeling the effect of the ion transport on fluid motion. With the multiphysics algorithm presented in this article,
the resulting multiply coupled, interacting system can be simulated efficiently on massively parallel supercomputers.
This algorithm is implemented in the waLBerla framework, whose modular software structure naturally supports
multiphysics simulations by allowing to flexibly combine different models. The largest simulation of the complete
system reported here performs more than 70 000 time steps on more than five billion (5 × 109 ) mesh cells for both
the hydrodynamics, as represented by a D3Q19 lattice Boltzmann automaton, and the scalar electric field. The
computations are executed in a fully scalable fashion on up to 8192 processor cores of a current supercomputer.
c 2017 Published by Elsevier Ltd.
Keywords: Parallel simulation; Electrokinetic flow; Electrophoresis; Fluid-particle interaction; MPI.

1. Introduction
1.1. Motivation
The motion of charged particles in fluids under the influence of electric fields occurs in a wide range
of industrial, medical, and biological processes. When the charged particles are immersed in liquids, their
migration caused by electric fields is termed electrophoresis. Due to the complex interplay of the physical
effects involved in such particle-laden electrokinetic flows, numerical simulations are required to analyze,
∗ Corresponding

author
Email address: dominik.bartuschat@cs.fau.de (Dominik Bartuschat)

1

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

2

predict, and optimize the behavior of these processes. To this end, we present a parallel multiphysics
algorithm for direct numerical simulations of electrophoretic particle motion.
Industrial applications that involve electrophoretic effects are electrofiltration [1, 2, 3] and electrodewatering [4]. Moreover, electrophoresis is utilized in electrophoretic deposition techniques for fabricating
advanced materials [5] and especially ceramic coatings [6, 7] in material science. Electrophoresis and electric
fields are also applied in many medical and biological applications. The trend towards miniaturization of
analysis processes has lead to the development of micro total analysis systems. Due to their high portability,
reduced costs, fast operation, and high sensitivity [8, 9], the design of such lab-on-a-chip systems has been a
highly active area of research for many years. These microfluidic systems require only small samples of liquid
and particles, which are transported, manipulated, and analyzed in structures of length scales from several
nm to 100 µm. Therefore, microfluidic separation and sorting of particles and cells are important steps of
diagnostics in such systems [10, 11]. Many of the employed techniques utilize electric fields to manipulate,
separate, and sort biological particles and macromolecules [8], such as cells [9, 10] or DNA [12].
At the small scales of microfluidic analysis systems, flow measurements are difficult or even impossible. Moreover, the complex coupling of hydrodynamic and electrostatic effects involved in electrophoretic
processes make predictions of electrophoretic motion challenging, especially for large numbers of particles.
Therefore, numerical simulations are essential to aid the design and optimization of electrophoretic systems.
The different physical effects in electrophoretic deposition can be better understood from insight gained
in simulations. By means of such simulations, electrophoretic sorting in lab-on-a-chip systems can be optimized for maximal throughput, sorting efficiency, and sorting resolution. A review of simulation methods for
electrophoretic separation of macromolecules is given in [13]. Also industrial applications of electrophoretic
deposition can be optimized with the help of simulations, as presented in [14] for a coating process.
1.2. Multiphysics Coupling Strategy
For simulations of electrokinetic flows with electrophoretic particle motion, the coupling between three
system components must be modeled: charged objects, fluid flow, and electric effects. The interacting
physical effects are sketched in Fig. 1. In the simulation method introduced in this article, the motion of the

ce
on
io

ty

ns

io
ot

e

si

m

rg

n
de

n
n

a
ch

r
fo

a
st

e
rc
fo

io

el

tro
ec

tic

Electroquasi-statics

object motion

Rigid body
dynamics

Fluid
dynamics
hydrodynamic force

Figure 1: Coupled physical effects of electrophoresis simulated with waLBerla and pe.

rigid, charged particles is modeled with Newtonian mechanics. The motion of the surrounding fluid, which
exerts hydrodynamic forces on the particles, is described by the incompressible Navier-Stokes equation. To
capture fluid-particle interactions, the hydrodynamic forces and the influence of the particle motion on the
fluid are modeled, based on the momentum exchange between fluid and particles. In this way, long-range
hydrodynamic interactions between individual particles and between particles and walls are recovered.
Moreover, electrostatic forces exerted by applied electric fields on the charged particles are modeled, which
cause the electrophoretic motion. The varying positions of the charged particles in return affect the electric
2

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

3

potential distribution in the simulation domain, based on their surface charge. Such a charge is carried by
most biomolecules such as cells, proteins, and DNA [15]. In fact, most substances acquire surface charges
when they get into contact with an aqueous medium [16] or electrolyte solution [17].
Electrostatic forces additionally act on the ions in the fluid and affect the fluid motion via body forces
present in locations of net charge. This net charge originates from the repulsion of co-ions in the fluid of
the same polarity as the surface charge and from the attraction of counter-ions. The ions in the fluid are
transported with the fluid flow, which in turn alters the electric potential distribution. In general, the motion
of ions in the fluid due to an electric field is governed by the Nernst-Planck equation, an advection-diffusion
equation that employs a continuum description and treats the ions as point charges. For flows in which
diffusion strongly dominates over advection and therefore quasi-thermodynamic equilibrium can be assumed
to hold, as considered in this paper, the electric potential due to the ion charge distribution is governed by
the Poisson-Boltzmann equation.
The region of the particle’s surface charge and of excess counter-ions in the fluid is denoted as electric
double layer (EDL). According to the Stern model [18] employed in this article, this double layer comprises
a region of ions attached to the surface and a diffuse part in which the ion concentration follows a Boltzmann distribution. At the surface of shear between the particle and the surrounding diffuse double layer,
the characteristic ζ-potential is defined. The employed equilibrium considerations based on the PoissonBoltzmann equation capture the dominant retardation effect of electrophoretic motion. This effect describes
the retardation of the charged particle motion by the action of the applied electric field on the opposite net
charge in the surrounding EDL. At high ζ-potentials, additionally the weaker relaxation effect occurs [16]
that is caused by a distortion of the EDL and can be captured by the Nernst-Planck equation.
In this article, we present an efficient parallel multiphysics simulation algorithm for electrophoresis on
a fixed Eulerian grid with a Lagrangian representation of moving particles. The particles are represented
by the physics engine pe [19, 20] as geometrically fully resolved three-dimensional objects. Dependent
on the electrostatic and hydrodynamic forces acting on the particles, the pe computes their trajectories
by rigid body dynamics and additionally resolves particle collisions. The pe is coupled to waLBerla
[21, 22, 23], a massively parallel simulation framework for fluid flow applications that employ the lattice
Boltzmann method (LBM) [24, 25]. By means of a modular software structure that avoids dependencies
between modules, functionality from different modules can be combined flexibly. To model fluid-particle
interactions, the LBM momentum exchange method [26, 27] implemented in waLBerla [28] is applied. For
the electrophoresis simulations, the LBM is performed with the two-relaxation time collision operator [29, 30]
and an appropriate forcing term for the electric body force due to the ions in the fluid. The electric potential
is represented by the Debye-Hückel approximation of the Poisson-Boltzmann equation that is discretized
with finite volumes whose mesh structure naturally conforms to the lattice Boltzmann (LB) grid. This
discretization also facilitates accommodating variable and discontinuous dielectricity values that vary in
time according to the particle positions, as required for simulating dielectrophoretic effects. By means of the
waLBerla solver module introduced in [31] together with the cell-centered multigrid solver implemented
therein, the resulting linear system of equations is solved. Since the counter-ions lead to a quicker decay of
the electric potential compared to the long-range electric potentials modeled in [31], the parallel successive
over-relaxation (SOR) method implemented in this module is an adequate choice.
In a previous article, we have shown that the implemented fluid-particle interaction algorithm for arbitrarily shaped objects efficiently recovers hydrodynamic interactions also for elongated particles [32]. Moreover,
we presented a parallel multiphysics algorithm for charged particles in the absence of ions in the fluid in
[31]. We have therein shown that several millions of charged particles with long-range hydrodynamic and
electrostatic interactions can be simulated with excellent parallel performance on a supercomputer. The
present paper extends these simulation algorithms by also considering ions in the EDL around the particles and their effect on the fluid motion, and presenting suitable parallel coupling techniques. Together
with the full four-way coupling of the fluid-particle interaction [31], the coupling with the quasi-equilibrium
representation of the electric effects results in a 7.5-way interaction.

3

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

4

1.3. Related Work
In the following, we give an overview of numerical methods for simulations of electrophoretic phenomena
that have been developed for different resolution levels. At the coarsest modeling scale, both fluid and
solid phase are described by an Eulerian approach. These continuum models represent the charged species
in terms of concentrations and are well suited for simulations with large numbers of unresolved particles.
In [33] two-dimensional electrophoresis simulations of biomolecule concentrations are presented. These
finite element simulations consider the effect of reactive surfaces and the electric double layer is included
through the biomolecule charge. Also three-dimensional parallel simulations of electrophoretic separation
with continuum approaches have been reported. The finite element simulations in [34] consider the buffer
composition and the ζ-potential at channel walls. In [35, 36] mixed finite element and finite difference
simulations of protein separation were performed on up to 32 processes.
At the finest level of resolution, fluid, ions, and particles are simulated by Lagrangian approaches.
These explicit solvent methods [37] typically apply coarse-grained molecular dynamics (MD) models to
describe the motion of fluid molecules and incorporate Brownian motion [38]. The mesoscale dissipative
particle dynamics method is applied in [39] to simulate electrophoresis of a polyelectrolyte in a nanochannel.
Another explicit solvent method is presented in [40] for simulating DNA electrophoresis, modeling DNA as
a polymer. In both methods the polymer is represented by bead-spring chains with beads represented by
a truncated Lennard-Jones potential and connected by elastic spring potentials. Explicit solvent models,
however, are computationally very expensive, especially for large numbers of fluid molecules due to pairwise
interactions [40]. Moreover, the resolution of solvent, macromolecules, and ions on the same scale limits the
maximal problem sizes that can be simulated [41]. Also the mapping of measurable properties from colloidal
suspensions to these particle-based methods is problematic [37].
The high computational effort is significantly reduced in implicit particle-based methods that incorporate
hydrodynamic interactions into the inter-particle forces. Such methods are applied in [37] and [42] to simulate
electrophoretic deposition under consideration of Brownian motion and van der Waals forces. Nevertheless,
these methods are restricted to few particle shapes and hydrodynamic interactions in Stokes flow.
Euler-Lagrange methods constitute the intermediate level of resolution. These approaches employ Eulerian methods to simulate the fluid phase, whereas the motion of individual particles is described by Newtonian mechanics. For simulations of particles in steady-state motion, the resolved particles can be modeled
as fixed while the moving fluid is modeled by an Eulerian approach. In [43] the finite volume method is
applied to simulate electrophoresis of up to two stagnant toroids in a moving fluid, employing the Hückel
approximation for a fixed ζ-potential and different electrical double layer thicknesses. The steady-state electrophoretic motion of particles with low surface potentials under a weak applied electric field in a charged
cylindrical pore is simulated in [44] for a single cylinder and in [45] for two identical spheres. In both cases,
a two-dimensional simulation with a finite element method for Stokes flow and Hückel approximation is
performed, exploiting the axial symmetry of the problem.
For electrophoresis at steady state perturbation approaches can be employed that are based on the assumption that the double layer is only slightly distorted from the equilibrium distribution for weak applied
electric fields (w.r.t. the field in the EDL). In addition to the equilibrium description based on the PoissonBoltzmann equation, small perturbations in the equilibrium EDL are considered in terms of linear correction
terms in the applied electric field for the ion distribution and the electric potential (see e. g. [46]). Using a
perturbation approach with finite elements, the electrophoresis of two identical spheres along the symmetry
axis of a cylindrical domain at pseudo steady-state were studied in [47]. Additionally to the hydrodynamic
and electric interactions, these axisymmetric simulations consider van der Waals forces for particles in close
proximity. In [48] a perturbation approach is applied to simulate a single colloid in a rest frame with periodic boundary conditions. The zeroth-order perturbation corresponds to the Poisson-Boltzmann equation,
which is solved by a constrained variational approach suggested in [49]. For the first-order perturbation,
the stationary Stokes equation is solved by a surface element method, the convection-diffusion for ionic
concentrations by a finite volume solver, and the Poisson equation by a fast Fourier Transform method.
More sophisticated Euler-Lagrange methods include direct numerical simulation (DNS) models that
represent the moving particles as geometrically fully resolved objects. These methods for particulate flows
4

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

5

include approaches with body-fitted moving meshes and fixed meshes. Moving meshes can be represented by
the Arbitrary-Lagrangian-Eulerian (ALE) formulation [50, 51] that employs moving unstructured meshes for
fluid-particle interaction problems. Such an ALE method is applied in [52] to simulate electrophoresis of a
single particle surrounded by a thin electrical double layer. The moving-mesh techniques require re-meshing
when the distortion of the updated mesh becomes too high, and subsequent projection of the solution onto
the new mesh. This overhead is circumvented in fixed-mesh techniques that allow the use of regular grids
and therefore the application of efficient solvers. The fluid particle dynamics method [53] falls into the latter
category, solving the Navier-Stokes and continuity equation on a fixed lattice, and representing the moving
solid particles by fluid particles of very high viscosity. By means of a concentration field that represents the
particle distribution, the particles affect the fluid viscosity and, together with forces acting on the particles,
the body force term of the Navier-Stokes equation. The rigid particles are geometrically modeled by the
Lennard-Jones potential and their motion is described by Newtonian mechanics [53]. This technique is
applied to simulate electrophoretic deposition of two charged particles in [54] and electrophoretic separation
in [55]. In these simulations the electrostatic interactions of the particles are modeled in terms of the body
force field, together with the advection and diffusion of ions and the resulting effect of the applied electric field
on the fluid motion. With this method, however, the particle rigidity is imposed by very high viscosity values
that restrict the time step size [56] and the Lennard-Jones potential restricts the particles shapes to spheres.
The smoothed profile method (SPM) [56] circumvents this time-step constraint by directly modeling the
particles as solid objects. Inside the particles and at solid-fluid boundaries that are represented by diffuse
interfaces, a body force is imposed on the fluid to model the effect of the particle motion on the fluid.
The fluid is again modeled on a fixed Cartesian grid and the particle motion with Newtonian mechanics,
where particle overlaps are typically prevented by a truncated Lennard-Jones potential. With this method,
electrophoresis of charged spherical particles is simulated in [57] for a constant, uniform electric field and
in [58] for an oscillating electric field. In both articles, the ion number concentration is modeled by an
advection-diffusion equation to recover non-equilibrium double layer effects. The SPM is also applied in [59]
to simulate electrophoresis of single cylinders and microtubules, employing the equilibrium representation of
the EDL. A further fixed-mesh technique is the immersed boundary method [60, 61, 62], where the rigid body
motion is imposed on the flow by body forces applied at the particle boundaries. This method, combined
with a finite volume method for solving the steady-state Poisson-Nernst-Planck equation system, is applied
in [63] to simulate the electrophoretic motion of up to three spherical particles in a two-dimensional setup.
Lattice-Boltzmann based methods are very well suited for parallel direct numerical simulations of fluidparticle interactions on fixed Cartesian grids. Both the Lagrangian particles and ions are often explicitly
modeled by molecular dynamics approaches that represent the rigid objects by repulsive potentials. In [64]
the electrophoresis of a colloidal sphere immersed in a fluid with counter-ions is simulated, modeling the
solvent by a lattice Boltzmann method. The charged sphere modeled with molecular dynamics is represented
by a raspberry model that comprised several beads connected by the finitely extensible nonlinear-elastic
(FENE) potential. Using a modified raspberry model with two spherical shells of beads solidly attached
to a larger spherical particle, this method is extended in [65] to simulate the electrophoresis of a spherical
Janus particle. The partially uncharged particle is surrounded by anions and cations represented by charged
beads. Electrophoresis simulations for a single highly charged spherical macro-ion in an electrolyte solution
with explicitly modeled positive and negative micro-ions are presented in [66]. Since the coupling of fluid
and macro-ion is performed via several particle boundary points, a single spherical particle is sufficient
to represent the macro-ion. A similar LB-MD method is applied in [67] to simulate the stretching of a
charged polyelectrolyte between parallel plates. The polyelectrolyte immersed in a liquid with explicitly
modeled counter- and co-ions is modeled by beads bonded together by the FENE potential. In all these
LB-MD simulations with explicitly modeled ions, hydrodynamic interactions are simulated with the LBM
and thermal fluctuations are added to both the fluid and the MD objects. The high computational effort
for modeling each individual ion by means of molecular dynamics, however, restricts the maximum feasible
problem size. With these approaches, only a limited number of ions per colloidal particle can be simulated
and the colloid radius is typically restricted to one order of magnitude larger than the ion size [48]. Therefore,
approaches based on continuum descriptions of the suspended ions are more practical for simulations of many
or for larger charged particles.
5

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

6

Alternatively to the continuum approach based on the Poisson-Boltzmann equation employed in this
article, electrophoresis can be simulated with the link-flux method [68] that models the advection and
diffusion of ions in terms of Nernst-Planck equations. The link-flux method employs the LBM for fluid
dynamics and models ion motion in terms of fluxes between lattice cells. In [69] this method is compared to
a LB-Poisson-Boltzmann approach for a fixed spherical particle in a periodic three-dimensional domain. The
aim is to examine the influence of particle motion and counter-ion concentration on the ζ-potential, leading to
the conclusion that for weakly perturbing electric fields or low Péclet numbers the equilibrium and dynamic
ζ-potentials are indistinguishable. The link-flux method is extended in [41] to support moving particles in
combination with the LB momentum exchange method. To ensure charge conservation in electrophoresis
simulations, appropriate moving boundary conditions for the solute fluxes are introduced. This method is
verified in [41] by electrophoresis simulations of up to eight particles.
1.4. Objectives and Outline
The primary goal of this paper is the introduction of a parallel multiphysics algorithm for electrophoresis
simulations together with validations of the physical correctness of the coupled algorithm for different particle
sizes. For this algorithm, waLBerla is augmented by an efficient boundary handling method that is able to
treat electric potential boundary conditions on the moving particles. Moreover, a joint parameterization for
the different coupled numerical methods is introduced. To achieve excellent computational performance, a
matrix-free representation of the linear system based on a stencil paradigm is used in the solver module [31].
For the linear Debye-Hückel approximation it is systematically exploited that these stencils are almost
uniformly identical throughout the simulation domain. The validation runs were performed on up to 8192
parallel processes of a modern supercomputer. Moreover, simulation results for the electrophoretic motion of
a single particle in a microchannel are presented, including visualizations of the electric potential distribution
and of the resulting flow field around the particle.
The equilibrium considerations in the present paper recover the predominant retardation effect due to
an opposing electrostatic force on the net opposite charge in the electrical double layer that counteracts the
particle motion. For the presented method, a computationally cheap and flexible SOR method is sufficient
to solve the electric potential equations. With our approach we aim for simulations of millions of charged
particles as in [31]. For these large numbers of particles the dynamics of an electrical double layer as in [41]
is computationally too expensive, even on modern supercomputers.
This paper is structured as follows: The physical background of fluid-particle interactions and electrophoresis are described in Sec. 2 and Sec. 3, respectively. In Sec. 4, the employed LB-momentum exchange
method for fluid-particle interactions is outlined, together with the finite volume discretization and the common parameterization concept for the coupled multiphysics methods. Then the extension of the waLBerla
framework for the electrophoresis algorithm is described in Sec. 5. Finally validation results for the electrophoretic motion of a spherical particle and visualizations of the resulting flow field and electric potential
distribution are presented in Sec. 6 before conclusions are drawn in Sec. 7.
2. Fluid-Particle Interaction
The macroscopic description of fluid behavior is based on the continuum hypothesis (cf. Batchelor [70])
that allows to consider a fluid as continuum, irrespective of the underlying molecular structure. In this case,
fluid properties can be represented by macroscopic quantities like density ρf , velocity ~u, and pressure p,
as functions of space and time. In terms of these quantities, fluid dynamics is described by conservation
laws for mass, momentum, and energy. In this article isothermal flows are considered, and therefore the
energy equation does not have to be solved. Moreover, non-continuum effects that become relevant for gas
flows at very small scales [71] are assumed to be negligible. Therefore, no-slip boundary conditions are
assumed to hold at solid-fluid interfaces, and slip velocities due to non-continuum Knudsen layer effects are
not considered.
Conservation of mass is described by the continuity equation. This equation can be derived by considering
a fixed control volume in the fluid relative to a stationary observer (Eulerian view). For incompressible fluids
6

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

7

the density of a fluid element is not affected by pressure changes [70, 72]. Then the density is spatially and
temporally constant (i. e. ρf = const), and the continuity equation reads as
∇· ~u = 0.

(1)

Conservation of momentum in a viscous, compressible fluid can be described in terms of the momentum
flux density tensor Π [72]. The temporal change of momentum in a control volume is balanced by the net
momentum flux through the surface of this volume and by external body forces f~b acting on the volume as
∂ (ρf ~u)
= −∇· Π + f~b .
∂t

(2)

The second-order tensor Π comprises a term for the convective transport of momentum and the total stress
tensor σ for the momentum transfer due to pressure and viscosity
Π = ρf ~u~u> − σ.

(3)

The total stress tensor σ can be decomposed into a part representing normal stresses related to the pressure
and a viscous part related to shear stresses. For incompressible Newtonian fluids, the stress tensor reads as


(4)
σ = −pI + µf ∇~u + (∇~u)> ,
where the first term with the second-rank identity tensor I contains the thermodynamic pressure p defined
according to Landau & Lifshitz [72] as used in the LBM literature [73]. The second term with dynamic
viscosity µf represents the shear stresses that are proportional to the rate of deformation [74] and result from
molecular transport of momentum [75]. With this stress tensor the incompressible Navier-Stokes equation
results from Eqn. (2), together with basic vector calculus and the continuity equation for compressible
fluids [76], as


∂~u
ρf
+ (~u · ∇)~u =
−∇p
+ µf ∆ ~u +
f~b
.
(5)
| {z }
|{z}
∂t
| {z }
|
{z
} pressure stress viscous stress external body force
inertial forces

This equation describes the balance of momentum change and the net force acting on a control volume
in terms of Newton’s second law. The left-hand side represents inertial forces acting on a fluid volume.
It comprises a term for local change of velocity and a term for convective acceleration i. e. the change of
velocity in space [77]. The right-hand side represents surface forces and body forces acting on the fluid
volume. Surface forces are short-range forces acting on the surface of the fluid element and are equivalent
to stress in the fluid [70]. Body forces, such as gravity or electrostatic force, act on the center of mass and
are represented by the force per unit volume of fluid element (or force density) f~b .
An important dimensionless quantity to characterize fluid flows is the Reynolds number Re = UνfL , with
µ
kinematic viscosity νf = ρff , characteristic velocity U , and characteristic length scale L. Flows in the regime
of creeping motion, where Re  1 and thus inertial forces are negligible, are termed Stokes flow.
The Stokes equations for incompressible Newtonian fluids resulting from Eqns. (5) and (1) read as
−∇p + µf ∆ ~u + f~b = 0,
∇· ~u = 0.

(6)

These equations are linear in both, velocity and pressure. Due to the linearity the superposition principle
holds, which is often utilized for fluid-particle interaction in Stokes flow and is employed for the validations
in Sec. 6.3.
Particles immersed in a fluid experience a force in case of relative fluid motion or a pressure gradient
in the fluid (e. g. due to gravity). This force exerted by the surrounding fluid can be calculated from the
7

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

8

stresses in the fluid next to the particle by integrating the stress tensor σ over the particle surface Γp [74, 46]
Z
F~part = σ · ~n dA.
(7)
Γp

Here, dA denotes the surface area elements and ~n the associated normal vector pointing into the fluid.
For simple cases, such as spherical bodies moving in unbounded Stokes flow, or single bodies moving in a
confined domain, analytical solutions for the particle motion are known as described in the following. More
complex fluid-structure interaction problems must be solved numerically, e. g. by the LBM.
For the computation of particle motion in incompressible fluids, hydrostatic effects that result e. g. from
gravitation acting on the fluid, do not have to be considered explicitly in the momentum equation (if the
buoyancy and gravitational force are directly applied to the particle). In this case, the force exerted by the
fluid on the particle according to Eqns. (7) and (4) is the drag force given by
Z
Z
~
~
~
Fd = − p dA + µf ∇~u + (∇~u)> dA,
Γp

Γp

where p is the hydrodynamic part of the total pressure.
The resistance to the motion of a sphere in an unbounded fluid at very low Reynolds numbers can be
calculated analytically from the above expression for the drag force and the Stokes equations (6) that govern
the fluid flow w.r.t. the imposed boundary conditions (BCs). For a sphere located at ~x = ~0, the unbounded
fluid is represented by the BC
~u → 0 as ~x → ∞
imposed on the fluid velocity in the Stokes equations (6). The resulting drag force acting on a rigid sphere
~ was derived by Stokes [78] as
of radius R that moves at constant velocity U
~.
F~d = −6πµf RU

(8)

This equation is commonly referred to as Stokes’ law.
For a sphere moving in a fluid subject to a constant force F~ , Stokes’ law relates the terminal steady-state
velocity of the sphere to the drag force exerted by the fluid. Such a constant force may e. g. be the Coulomb
force that acts on a charged particle in an electric field. The terminal sphere velocity is then obtained from
the balance of the external force and the drag force F~ + F~d = ~0 as
~ =
U

1
F~ .
6πµf R

(9)

A particle moving in a confined domain experiences a retardation caused by surrounding walls. Consequently, the drag force on a sphere that moves in a viscous fluid limited by walls is higher than the force
according to Stokes’ law. The effect of walls on a moving particle in Stokes regime can be determined by
means of the method of reflections, as described in detail in [74]. Happel & Bart [79] employed this method
to obtain a first-order correction to the drag force on a sphere settling in a long square duct with no-slip
walls. Miyamura et al. [80] found polynomial expressions for the increased drag by fitting the coefficients to
experimentally obtained settling velocities of spheres in different confining geometries. The correctness of
the wall effect recovered in LBM simulations with the fluid-particle interaction algorithm employed in this
article was verified in [76] against these expressions.
3. Electrokinetic Flows
The transport of ions in fluids subject to electric fields that occurs in electrokinetic flows can be modeled
by means of a continuum theory, similar to the description of fluid dynamics by the Navier-Stokes equation.
8

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

9

Instead of modeling individual ions and their interactions, local ion concentrations ni and fluxes ~ji of the
different ionic species i are considered. Based on these macroscopic quantities the ion transport in dilute
electrolyte can be described by the Nernst-Planck equation combined with the law for the conservation of
ionic species in a solution
∂ni
= −∇· ~ji
(10)
∂t
in the absence of chemical reactions. Here ni denotes the number density or concentration that is related
to the molar concentration as ci = NnAi with Avogadro’s number NA . The total ionic flux ~ji of species i
comprises an advective flux with a common mass average velocity ~u for all species and fluxes relative to the
advective flux due to diffusion and electric migration [46]. This relation is expressed by the Nernst-Planck
equation
~ji =
ni ~u
− Di ∇ni − ni µ∗i ∇Φ,
(11)
|{z}
| {z }
| {z }
advective flux

diffusive flux

migration flux

where Di and
represent the spatially homogeneous diffusion coefficient and ionic mobility of species i,
respectively, and ∇Φ the local electric potential gradient. The ionic mobility is defined as µ∗i =: DkiBzTi e ,
where e denotes the elementary charge, zi the valence of a given ionic species, kB the Boltzmann constant
and T the temperature.
To model the influence of the charged ions on the electric potential governed by the Poisson equation
µ∗i

− ∆ Φ(~x) =

ρe
εe

(12)

for spatially uniform fluid permittivity εe , the ion charge distribution is considered in terms of the local
mean macroscopic charge density as
X
ρe =
e zi ni .
(13)
i

The Poisson-Nernst-Planck equation system Eqns. (10)–(12) is highly nonlinear, and solving the overall system is computationally very expensive, especially for electrophoresis of many particles. Therefore
the problem is simplified by restriction to equilibrium considerations based on the Boltzmann distribution
that capture the dominant electrophoretic effects. The resulting Poisson-Boltzmann equation holds for
(quasi-)thermodynamic equilibrium when the ion distribution is not affected by fluid flow or by externally
applied electric fields. Therefore the electric potential ψ resulting from the non-uniform ion distribution in
the EDL is considered in the following, instead of the total electric potential Φ = ψ + ϕ that additionally
comprises the potential ϕ of the externally applied electric field.
The Boltzmann distribution for ions can be derived from the Nernst-Planck equation, as outlined in [46].
Considering the Nernst-Planck equation (11) in one dimension and at equilibrium, i. e., for zero macroscopic
fluid velocity u = 0 and ionic flux ji = 0, results in
d ni
zi e d ψ
= −ni
dx
kB T dx

(14)

for the above definition of µ∗i . In this case, the Péclet number P e = UDL for mass transfer relating the
advection rate to diffusion rate becomes zero. Here U is the fluid speed, L a characteristic length scale, and
D the diffusion coefficient. Applying the chain rule to the left-hand side of Eqn. (14) and integrating from
a reference point in the bulk with potential ψ∞ and concentration ni∞ , yields
zi e(ψ − ψ∞ )
kB T
ni = ni∞ e
.
−

(15)

Setting the reference potential ψ∞ in the electroneutral bulk solution to zero recovers the Boltzmann distribution with the number density ni∞ at the location of the neutral state.
9

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

10

From Poisson’s equation (12) with net charge density according to Eqn. (13) and the obtained Boltzmann
distribution, the Poisson-Boltzmann equation follows as
zi e ψ
−
e X
− ∆ψ =
zi ni∞ e kB T ,
εe i

(16)

relating the electric potential ψ to the ion concentrations at equilibrium. For binary, symmetric electrolyte
solutions comprising two species of valence z = −z− = z+ , the Poisson-Boltzmann equation takes the form


2 z e n∞
zeψ
− ∆ψ = −
sinh
.
(17)
εe
kB T
For low ζ-potentials compared to the thermal voltage kB T /e, the term kzBe ψT in Eqn. (17) becomes
smaller than unity. At room temperature this is fulfilled for ζ  25.7zmV [81]. In this case the approximation
sinh(x) ≈ x is accurate, up to a small error of order O(x3 ) (cf. Taylor’s expansion). With this linearization,
the symmetric Poisson-Boltzmann equation simplifies to the Debye-Hückel approximation (DHA)
− ∆ψ = −

2 e2 z 2 n∞
ψ = −κ2 ψ.
εr ε0 kB T

(18)

This equation was originally derived by Debye & Hückel [82] for strong electrolytes [81]. The parameter κ,
defined by
r
εr ε0 kB T
κ :=
,
(19)
2 e2 z 2 n∞
is commonly referred to as Debye-Hückel parameter. Moreover, the charge density in the fluid is then given
by
ρe (ψ) = −κ2 εe ψ.
(20)
In this article, we consider spherical particles with uniform ζ-potential distribution as depicted in Fig. 2.
The electric potential ψ for such a particle of radius R is represented by the Debye-Hückel equation in

E∞
r
ψ =ζ

~u = −U

θ

R

x

λD
~ex
Figure 2: Electrophoresis setup of a stationary (negatively) charged sphere of radius R, surrounded by a double layer and
subject to an applied electric field in opposing fluid flow. Similar to [16].

spherical-polar coordinates as

1 d
r2 dr



r2

dψ
dr
10



= κ2 ψ,

(21)

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

11

with radial distance r from the sphere center and subject to the Dirichlet BCs
ψ=ζ
ψ→0

at
as

r = R,
r → ∞.

(22)

Solving this equation subject to these BCs results in the electric potential distribution in the surrounding
EDL (and beyond) as [81]
R
ψ(r) = ζ e−κ(r−R)
for r ≥ R.
(23)
r
For the Debye-Hückel approximation, an analytical solution for the electrophoretic velocity of a single
spherical particle at steady state and arbitrary EDL thickness has been derived by Henry [83]. Initially the
problem was formulated to account for finite conductivities of particle and medium, and the potential in
the EDL was described by the Poisson-Boltzmann equation. The final results, however, were provided for
insulating spheres, and sufficiently low ζ-potentials for the Debye-Hückel approximation to hold. Therefore
the Debye-Hückel approximation is employed in the following to represent the ion distribution around the
particles under electrophoretic motion.
Instead of modeling a sphere moving at terminal speed U , Henry [83] considered the analogous problem of
a stationary sphere in a steadily moving liquid. To this end, the opposing velocity −U was imposed on the
overall system by setting ~u r→∞ = −U ~ex in the liquid far from the particle. As shown in Fig. 2, a sphericalpolar coordinate system fixed at the particle center with radial distance r and polar angle θ was used. Under
the assumption that the electric potential in the EDL is not distorted from its equilibrium distribution by
the applied field and the fluid flow, the potentials ϕ and ψ were linearly superimposed. Therefore, the
electric potential ψ in the diffuse double layer is described by the Poisson-Boltzmann equation and the
applied potential ϕ by a Laplace equation. The BCs for the Laplace equation applied by Henry represent
the insulating particle by homogeneous Neumann BCs at the particle surface and impose the applied field
by the inhomogeneous Neumann condition ∂ϕ/∂x r→∞ = −E∞ . For the Poisson-Boltzmann equation, the
ζ-potential at the hydrodynamic radius R and the decaying potential were imposed as given in Eqn. (22).
Making use of the equations for the electric potential, the Stokes equations for steady-state creeping flow
with body force term on the right-hand side
−µf ∆ u + ∇p

= −ρe ∇ (ϕ + ψ)

∇· ~u =

(24)

0,

(25)

were solved by Henry [83]. In addition to the BC imposing the opposing velocity far from the particle to
bring the whole system to rest, the no-slip condition ~u r=R = 0 was applied at the particle surface.
From the flow field around the particle, the force acting on the particle was obtained by integrating
the normal stresses over the sphere surface. To the resulting force that comprises Stokes drag and electric
components, the electrostatic force on the particle due to its fixed surface charge was added. The total force
must vanish at steady motion and was thus equated with zero, resulting in the electrophoretic velocity



ZR
ZR
ψ
ψ
ε
~ ext
~ EP = e ψR + R3 5R2
dr − 2
dr E
(26)
U
µf
r6
r4
∞
∞
|
{z
}
=ζ f (κR)

obtained by Henry for an insulating particle. The function f (κR) introduced in [83] is usually referred to
as Henry’s function. In [84] the following expression is derived that approximates the integral equations as


f (κR) =

2

1 +
3

1

2 1+
11

2.5
κR (1 + 2e−κR )



3  ,


(27)

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

12

with a relative error below 1 % for all values of κR [84].
In terms of Henry’s function f (κR) [83], the electrophoretic mobility of a spherical, non-conducting
particle reads as
εe
ζ f (κR) .
(28)
µEP =
µf
This solution is correct to the first order of the ζ-potential, since the relaxation effect is neglected [85]. With
the definition of the electrophoretic mobility µEP := UEP /Eext as electrophoretic speed per unit applied field
the electrophoretic velocity of the particle can be obtained.
Henry’s analytical solution for the electrophoretic velocity of a spherical particle of radius R in an
unbounded electrolyte solution of dynamic viscosity µf with Debye-Hückel parameter κ, subject to an
~ ext reads as
applied field of strength E


~ EP = 2εe ζ 
U
1 +
3µf

1
h

2 1+

2.5
κR(1+2e−κR )

 ~
i3  E
ext ,

(29)

according to the expression for the electrophoretic mobility of Ohshima given in Eqn. (27).
For the electrophoresis simulations in this article, the particle charge must be known to compute the
electrostatic force on the particle. Analytical solutions for electrophoretic motion such as Henry’s equation,
however, are typically given in terms of the ζ-potential, which is defined at the slip surface between the
compact and diffuse EDL layer. Since the particle charge is acquired as a surface charge, for a given ζpotential the surface charge (density) enclosed by the slip surface is therefore needed in the simulations.
The surface charge density is hereby obtained from the overall surface charge bound at the fluid-particle
interface and in the Stern layer. This approach is justified by the fact that the electric potential at the Stern
surface and the ζ-potential can in general be assumed to be identical [86].
The relation of the surface density σs to the ζ-potential is obtained from the Neumann BC on the surface
of the insulating particle [85]
dψ
σs = −εe
(30)
d~r r=R
in case these electrical properties do not vary in angular direction. This condition holds for insignificant
permittivity of the insulating particle compared to the fluid permittivity εe . Alternatively, this relation can
be derived from the electroneutrality condition [46]. With the spatial distribution of ψ around the spherical
particle according to Eqn. (23) the ζ–σs relationship follows from Eqn. (30) as


qs
1 + κR
σs =
=
ζ
ε
.
(31)
e
4πR2
R
For a spherical particle with an EDL potential ψ described by the spherical symmetric Poisson-Boltzmann
equation, the more general ζ–σs relationship
v



 
u
zeζ
u
8 ln cosh

 zeζ u
2
1
2εe κkB T
1
4kB T
u


+



u1 +
sinh 
(32)
σs =

2
zeζ
zeζ
ze
κR
 2kB T t

(κR)
2
2
cosh
sinh
4kB T
2kB T
for 1-1 and 2-1 electrolyte solutions is derived in [87]. The relative error of this approximation w.r.t. the
exact numerical results computed by [88] is below 1 % for 0.5 ≤ κR < ∞ [85]. This relationship is applied in
the electrophoresis simulation validation in Sec. 6.3 to compute the particle charge for a given ζ-potential.
The applied ζ–σs relationship is derived for electric potentials governed by the spherical Poisson-Boltzmann
equation and is thus more general than the ζ–σs relationship (31) for the Debye-Hückel approximation. For
the simulation parameters used in Sec. 6.2, the deviation of σs for the general relationship w.r.t. the exact
value of the Debye-Hückel approximation is about 0.2 % and hence negligible.
12

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

13

4. Numerical Modeling
4.1. Lattice Boltzmann Method with Forcing
The LBM is a mesoscopic method for the numerical simulation of fluid dynamics based on kinetic theory
of gases. This method statistically describes the dynamics of ensembles of fluid molecules in terms of particle
distribution functions (PDFs) that represent the spatial and velocity distribution of molecules in phase space
over time. The temporal and spatial variation of PDFs, balanced by molecular collisions, is described by
the Boltzmann equation. The solution of this equation converges towards the Maxwell-Boltzmann velocity
distribution of molecules in local thermodynamic equilibrium. For small deviations from this equilibrium,
the Navier-Stokes equations can be derived from the Boltzmann equation by means of a Chapman-Enskog
expansion (cf. [89, 90]).
In the LBM, the phase space is discretized into a Cartesian lattice Ωdx ⊂ RD of dimension D with
spacing dx, and a set of Q discrete velocities ~cq ∈ RD , q ∈ {1, . . . , Q}. Moreover, time is discretized as
cq are chosen
Tdt = {tn : n = 0, 1, 2, . . .} ⊂ R+
0 , with a time increment of dt = tn+1 − tn . The velocities ~
such that within a time increment, molecules can move to adjacent lattice sites or stay at a site. Associated
with each of these velocities is a PDF fq : Ωdx × Tdt 7→ R. A forward-difference discretization in time and
an upwind discretization in space [91] result in the discrete lattice Boltzmann equation
fq (~xi + ~cq dt, tn + dt) − fq (~xi , tn ) = dtCq + dtFq ,

(33)

with lattice site positions ~xi , discrete collision operator dtCq , and discrete body-force term Fq . This equation
describes the advection of PDFs between neighboring lattice sites and subsequent collisions.
In general, the collision operator can be represented in terms of the collision matrix S as dtCq =
P
eq 
(cf. [92]), with the vector f~ := (f0 , f1 , . . . , fQ−1 )> of the PDFs fq , and f~eq of the equij Sqj fj − fj
librium distributions fqeq . The latter are obtained from a low Mach number expansion of the MaxwellBoltzmann distribution [89]. With a representation of the macroscopic fluid density as ρf = ρ0 + δρ in
terms of a reference density ρ0 and a density fluctuation δρ, the equilibrium distribution function for the
incompressible LBM is derived in [93] as



(~u · ~u)
(~cq · ~u) (~cq · ~u)2
+
−
,
(34)
fqeq (ρ0 , ~u) = wq ρf + ρ0
c2s
2c4s
2c2s
where ‘·’ denotes the standard Euclidean scalar product. This distribution function depends on ρf , the
macroscopic fluid velocity ~u, and lattice-model dependent weights wq . At each instant of time, ρf and ~u are
given by moments of PDFs as
P
ρf (~xi , t) = fq (~xi , t),
qP
(35)
~cq fq (~xi , t).
~u(~xi , t) = ρ10
q

c2s ρf (~xi , t)

according to the equation of √
state for an ideal
Moreover, the pressure p is given as p(~xi , t) =
gas. We employ the D3Q19 model of [94] with thermodynamic speed of sound cs = c/ 3 for the lattice
velocity c = dx/dt. For this model, the weights wq are: w1 = 1/3, w2,...,7 = 1/18, and w8,...,19 = 1/36. As
discussed in [93], fqeq recovers the incompressible Navier-Stokes equation with at least second-order accuracy
u|
2
of the Mach number Ma := |~
cs , O(Ma ). In LBM simulations, the kinematic fluid viscosity νf is generally
determined by a dimensionless relaxation time τ of the collision operator as


1 2
c dt.
(36)
ν= τ−
2 s
As shown in [95], with this definition the LBM is second-order accurate in space and time.
Among the different collision operators available for the LBM, we employ the two-relaxation-time (TRT)
collision operator of [29, 30]
X



Sqj fj − fjeq = λe fqe − fqeq,e + λo fqo − fqeq,o ,
(37)
j

13

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

14

with the relaxation parameters λe and λo for even- and odd-order non-conserved moments, respectively.
Alternative collision operators either have disadvantages regarding stability and accuracy [96] such as the
BGK model [97] or are computationally more costly such as the MRT operator [92] or the cumulant operator
[98]. To ensure stability, both relaxation parameters should be within the interval ] − 2, 0[, cf. [99, 30].
The even relaxation parameter is related to the dimensionless relaxation time by λe = − τ1 and therefore
determines the fluid viscosity. The free parameter λo is set to λo = −8(2 − 1/τ )/(8 − 1/τ ) in this article,
which prevents τ -dependent boundary locations for no-slip BCs as they arise for the BGK operator. Instead,
walls aligned with the lattice dimensions are fixed half-way between two lattice sites, as shown in [100]. For
the TRT operator, the PDFs are decomposed as fq = fqe + fqo into their even and odd components
fqe = 12 (fq + fq̄ )
fqo = 12 (fq − fq̄ )

and
and

fqeq,e = 12 (fqeq + fq̄eq )
fqeq,o = 12 (fqeq − fq̄eq ),

with ~cq = −~cq̄ . The local equilibrium distribution function in Eqn. (34) is then given by


ρ0
ρ0
u · ~u) + 2c
cq · ~u)2
fqeq,e = wq ρf − 2c
2 (~
4 (~
s
s
fqeq,o = wq ρc20 (~cq · ~u).

(38)

(39)

s

At each time step tn ∈ Tdt the lattice Boltzmann method performs a collide– and a stream step
f˜q (~xi , tn ) = fq (~xi , tn ) +λe [fqe (~xi , tn ) − fqeq,e (~xi , tn )]
+λo [fqo (~xi , tn ) − fqeq,o (~xi , tn )]

(40)

fq (~xi + ~eq , tn + dt) = f˜q (~xi , tn ) + dtFq ,

(41)

where f˜q denotes the post-collision state and ~eq = ~cq dt a discrete lattice direction.
In the stream step, the product of dt and the forcing term Fq is added to the post-collision PDFs. The
term Fq considers the external effect of body forces acting on the fluid. In this article, the discrete forcing
term according to Luo [101] is employed as


(~cq − ~u) (~cq · ~u)
~
+
~
c
Fq = w q
(42)
q · fb ,
c2s
c4s
with f~b representing the electrical body force per unit volume. The forcing terms lead to an additional
term in the continuity equation that arises for spatially varying external forces, as shown in [102, 103]. This
additional term can be removed by incorporating the external force in the momentum density definition as
!
1 X
dt ~
~u =
fq~cq + fb .
(43)
ρ0
2
q
Thus, [102] suggest to use the forcing term (42) in combination with the modified momentum density
definition in Eqn. (43) for the BGK. We therefore use this forcing term with second-order accuracy, together
with the modified momentum density for the resulting velocity ~u.
To increase the computational efficiency of the implementation, the compute-intensive collide step and
the memory-intensive stream step are fused to a stream-collide step. In the simulations presented in this
article, no-slip and free-slip BCs are applied as described in [32].
4.2. Momentum Exchange Approach
To model the fluid-particle interaction and the resulting hydrodynamic interactions of the particles, the
momentum exchange approach suggested by [27] is employed in this article. This approach computes the
momentum exchange between the fluid and the suspended rigid particles from PDFs surrounding these
solid objects. The implementation of this approach in waLBerla has recently been applied to simulate
fluid-particle interactions also at Reynolds numbers beyond the Stokes regime, as presented in [104, 105, 106].
14

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

15

For the momentum exchange approach, the solid particles are mapped onto the lattice by considering
each cell whose center is overlapped by a particles a moving obstacle cell. All other lattice cells are fluid cells,
resulting in a staircase approximation of the particle surfaces. These surfaces are represented by surface cells
indicated by subscript s in the following. On the fluid cells denoted by subscript F , the LBM is applied. To
model the momentum transfer from the particles to the fluid, the velocity bounce-back BC
ωq
fq̄ (~xF , tn + dt) = f˜q (~xF , tn ) − 2 2 ρ0~cq · ~us
cs

(44)

is applied at fluid cells with position ~xF = ~xs + ~eq̄ adjacent to a surface cell at ~xs . This boundary condition
introduced in [26] matches the fluid velocity to the local particle surface velocity ~us .
From the sum of all force contributions due to the momentum transfer from fluid cells to neighboring
surface cells (cf. [32]), the overall hydrodynamic force on the particle can be obtained according to [27] as

X X 
dx3
ωq
F~h =
2f˜q (~xF , tn ) − 2 2 ρ0~cq · ~us ~cq
.
(45)
cs
dt
s
q∈Ds

Here, Ds is the set of direction indices q in which a given particle surface cell s is accessed from an adjacent
~ h is given by substituting ~cq × (~xs − ~xC ) for the last ~cq term in
fluid cell. Analogously, the overall torque M
Eqn. (45), with ~xC representing the particle’s center of mass.
The mapping of the solid particles onto the lattice results in fluid cells appearing and disappearing as
the particles move. Therefore, at uncovered lattice sites the PDFs must be reconstructed. We set the PDFs
at those fluid cells to the equilibrium distribution f eq (ρ0 , ~us (~xs (tn − dt)) according to Eqn. (34) dependent
on the local particle surface velocity from the previous time step.
4.3. Finite Volume Discretization for Electric Potential Equations
To solve the Debye-Hückel approximation a cell-centered finite volume scheme is applied on the Cartesian
lattice Ωdx introduced in Sec. 4.1. Associated with this lattice of spacing dx that represents the computational
domain Ω ⊂ R3 is the three-dimensional cell-centered grid Gdx defined (cf. [107]) as




j − 1/2

V
V 
Gdx := ~xi ∈ Ω ~xi =  k − 1/2  dx i ∈ i , (j, k, m) ∈ J
,
(46)


m − 1/2
with Ωdx = Ω ∩ Gdx . For
V indexing of lattice cells by tuples (j, k, m) of cell indices in the three spatial
dimensions, the index set J := {(j, k, m) | j = 1, . . . , lx ; k = 1, . . . , ly ; m = 1, . . . , lz } is introduced.
Here
V
lx , ly , and lz represent the numbers
of
cells
in
x,
y,
and
z-direction,
respectively.
The
set
is
related
J
V
to the
V setVof single cell indices i := {i | i = 1, . . . , lx · ly · lz } used for the LBM by a bijective mapping
g : i → J , according to Eqn. (46).
The finite volume discretization of the Debye-Hückel approximation Eqn. (18) includes volume integration
over each lattice cell Ωi = Ωklm := [~xj−1,k,m , ~xj,k,m ] × [~xj,k−1,m , ~xj,k,m ] × [~xj,k,m−1 , ~xj,k,m ] and applying
the divergence theorem to the Laplace operator ∆ = ∇· ∇, resulting in
I
Z
2
~
−
∇ψ(~x) dΓi + κ
ψ(~x) d~x = 0 ∀ Ωi ∈ Ωdx .
(47)
∂Ωi

Ωi

Here, ∂Ωi denotes the closed surface of the cell, and ~Γi is a surface directed element. The cell surface
consists of six planar faces with constant outward unit normal vectors ~niq (q = 1, . . . , 6). Therefore, the
surface integral can be decomposed into a sum of integrals [108] as
I
−
∂Ωi

∇ψ(~x) d~Γi = −

Z
6
X
q=1 ∂Ω
iq

15

∇ψ(~x) · ~niq dΓiq ,

(48)

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

16

where q represents LBM direction indices introduced in Sec. 4.1, and ∂Ωiq is the common face with the neighboring cell in direction q. The gradients ∇ψ(~xi ) · ~niq in normal direction of the faces ∂Ωiq are approximated
at the face centers by central differences of ψ from a neighboring and the current cell as
∇ψ(~xi ) · ~niq

~
xi + 12 ~
eq

≈

ψ(~xi + ~eq ) − ψ(~xi )
.
dx

(49)

Here, ~eq represents the corresponding lattice direction introduced in Sec. 4.1.
Substituting the approximation of the normal fluxes across the surfaces of area dΓiq = dx2 , Eqn. (49) into
Eqn. (48) results in
I
6
X
ψ(~xi + ~eq ) − ψ(~xi ) 2
−
dx .
(50)
∇ψ(~x) d~Γi ≈ −
dx
q=1
∂Ωi

Applying the above finite volume discretization to the linear term of the Debye-Hückel approximation
results in
Z
κ2
ψ(~x) d~x ≈ κ2 ψ(~xi ) dx3 .
(51)
Ωi
DHA
This additional term enters the central element of the resulting seven-point stencil Ξdx
as




0
−1
0
0
0
0
0
0
0
1 
 −1 6 + κ2 dx2 −1 0 −1 0 ,
0
−1
0
dx2
0
0
0
0
0
0
0
−1
0

(52)

and the right-hand side for each unknown is zero.
4.4. Parametrization for Electrokinetic Flows
Numerical simulations are typically performed in terms of dimensionless parameters. To ensure consistently good numerical accuracy independent of the simulated scales, the physical quantities are mapped to a
computationally reasonable numerical value range. In LBM simulations usually the quantities are expressed
in lattice units. Therefore, physical values must be converted to lattice values before the simulation and
vice versa to obtain physical values from simulation results. In the following the lattice unit (LU) system
employed in this article is presented, providing a common parameterization also for further (electrokinetic)
flow scenarios [76].
Physical simulation parameters are given in terms of the international system of quantities (ISQ) associated to the international system of units (Système International d’Unités, SI [109]) [110]. The value of a
quantity is defined as product of a numerical value and a unit, e. g. 10−6 m. The SI system comprises the
base units displayed in Tab. 1, together with the corresponding mutually independent base quantities length,
time, mass, electric current, thermodynamic temperature, amount of substance, and luminous intensity.
Moreover, derived units such as N = kgs2m (Newton) are defined in the SI system as products of powers of
base units [109].
For LBM simulations the base quantities are length, time, and mass density. The corresponding lattice
Table 1: Physical and LBM base quantities for electrokinetic simulations.

kind of quantity

length

time

mass

electr.

thermodyn.

chem.

photometr.

ISQ quantity
SI unit

x
m

t
s

m
kg

I
A

T
K

n
mol

Iv
cd

LBM quantity
lattice unit

x
dx

t
dt

ρ
ρ0

Φ
V

–
–

–
–

–
–

16

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

17

base units are the (physical) lattice spacing dx, time increment dt, and fluid reference density ρ0 . Therefore,
the numerical values of these quantities in lattice units become unity, as shown in Tab. 2. The lattice
parameters representing these dimensionless numerical values in LUs are indicated with subscript L in
the following, e. g. dxL for the lattice spacing. Performing LBM computations on such normalized lattice
parameters saves numerical operations: Additional scaling factors are avoided, e. g. in the LBM streamcollide step, when computing the equilibrium distribution function given in Eqn. (34) (cf. cL = dxL /dtL = 1)
or the macroscopic velocity according to Eqn. (35) (cf. ρ0,L = 1).
For the electrokinetic simulations, the electric potential Φ is chosen as electric base quantity. This quantity, however, is not scaled on the lattice but keeps its numerical value that typically lies in a range not too
far from unity. In contrast to the SI system, the LU system for electrokinetic simulations requires no base
units corresponding to the thermodynamic temperature T or the amount of substance n:
In the simulations, temperature appears only in combination with the Boltzmann constant as energy,
i. e., E = kB T , with the derived unit [E] = [m] [x]2 /[t]2 . With the relation of mass and mass density
[m] = [ρ0 ] [x]3 , this unit can be represented in lattice base units (see Tab. 1).
Moreover, by representing the molar concentration with unit [ci ] = mol
l in terms of the number density as
1
ni = ci NA with Avogadro’s number NA = 6.022 14 · 1023 mol
, the unit ‘mol’ cancels out, yielding [ni ] = [x]−3 .
With the choice of the potential Φ as base quantity, the electric current I becomes a derived quantity
in the LU system. The unit of I can be derived from the energy in electrical units [E] = [I] [Φ] [t] and the
above definition of energy in terms of lattice units [E] = [ρ0 ] [x]5 /[t]2 . Equating both relations results in the
derived lattice unit
[ρ0 ] [x]5
.
(53)
[I] =
[Φ] [t]3
In Tab. 2 different physical quantities and their SI units are displayed, as well as their representation by
(dimensionless) lattice parameters and their lattice units. The conversion of physical quantities to lattice
Table 2: Relation of physical quantities and lattice parameters for electrokinetic simulations.

physical quantity

SI unit

dx (lattice spacing)
dt (time increment)

m
s

ρ0 (fluid density)
Φ (electr. potential)

kg
m3

L (length)
ν (kinem. viscosity)
~u (velocity)
m (mass)
F~ (force)
I (electr. current)

V

m

m2
s
m
s

kg

kg m
s2

A

lattice parameter
dxL
dtL
ρ0,L
ΦL
LL
νL
~uL
mL
F~L
IL

e (elem. charge)

As

eL

ε0 (vacuum permittivity)
~ (electr. field )
E

As
Vm
V
m

ε0,L
~L
E

J

EL

E (energy)

numerical value
1
dx dx (= 1)
1
dt dt (= 1)
1
ρ0 ρ0 (= 1)
1
V Φ
1
dx L
dt
ν
dx2
dt
u
dx ~
1
m
ρ0 dx3
dt2
F~
ρ0 dx4
3
V dt
I
ρ0 dx5
V dt2
e
ρ0 dx5
V dt2
ε
ρ0 dx5 0
dx ~
V E
dt2
E
ρ0 dx5

lattice unit
dx
dt
ρ0
V
dx

dx2
dt
dx
dt

ρ0 dx3
ρ0 dx4
dt2
ρ0 dx5
dt3 V
ρ0 dx5
dt2 V
ρ0 dx5
V dt2
V
dx
ρ0 dx5
dt2

units requires their division by powers of the LBM base quantities with the corresponding SI units. Since the
potential Φ has the same numerical value in physical and lattice units, the LBM base unit of the potential is
dx5
simply ‘1 V’. Therefore, the derived lattice unit ‘Ampere’ for the electric current is given by A = ρV0 dt
3 (see
Eqn. (53)). The corresponding scale factors for converting the physical parameters (e. g. ν) to the associated
17

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

18

lattice parameters (e. g. νL ) are shown in Tab. 2. Multiplication with the inverse scale factors converts the
lattice parameters back to physical quantities.
5. Extension of waLBerla for Electrophoresis Simulations
Electrophoresis simulations require the mutual coupling of fluid dynamics, rigid body dynamics, and
electro-statics, as shown in Fig. 1. In addition to electrostatic and hydrodynamic interactions, the applied
field acts on the EDL charge and thereby affects fluid flow and particle motion. For the electrophoresis
algorithm presented below, the equilibrium description of the EDL potential in terms of the linear DebyeHückel equation is employed. Therefore, the predominant retardation effect is recovered in the simulations.
The applied potential ϕ and the EDL potential ψ are linearly superimposed (cf. Henry’s equation, Sec. 3),
which is valid for weak applied fields when the EDL distortion by the field is negligible. Thus, the applied
electric field can be imposed directly, without solving the associated Laplace equation. In the following the
main concepts of the waLBerla framework are described, together with the functionality for electrophoresis
simulations implemented therein.
5.1. Design Concepts of waLBerla
WaLBerla is a framework for massively parallel multiphysics simulations with a MPI-based distributed
memory parallelization that is specifically designed for supercomputers. The main software design goals of
this framework are flexibility to combine models of different effects, extensibility to allow the incorporation
of further effects and details, e. g., for electrokinetic flows, and generality to support further applications [76].
These goals are reached by integrating the coupled simulation models into waLBerla in a modular fashion
that avoids unnecessary dependencies between the modules. This way, the modules can be augmented
by more sophisticated models or models tailored to a certain application, and functionality from different
modules can be combined flexibly. The modular code structure also provides excellent maintainability, since
modifications of the code in one module do not affect other modules. Developers can therefore efficiently
locate faulty modules and find bugs inside these modules, also by systematically utilizing automatic tests.
In addition to modules, the waLBerla code structure comprises a core for sequence control that initializes data structures, performs the time stepping, and finalizes the simulation on each parallel process.
By means of applications, multiphysics simulations can be defined by assembling the associated functionality and coupled models from the modules. The coupling strategy for multiphysics simulations is based
on accessing mutually dependent data structures (see [31] for more details). These data strucutres are
defined in the modules that implement models for the different physical effects. Also infrastructural and
utility functionality is encapsulated in modules, e. g., for domain setup, MPI communication, BC handling,
parameterization, or simulation data output. For parallel simulations the discretized simulation domain is
decomposed into equally sized blocks of cells that can be assigned to different parallel processes. In this block
concept, each block contains a layer of surrounding ghost cells that is needed for BC treatment and parallelization. For parallelization, cell data of neighboring processes is copied to the ghost layer, dependent on
the data dependencies of the unknowns located on a given block. Moreover, metadata of a block specifies its
location in the simulation domain or its rank for MPI communication. The communication concept provides
a simple and flexible communication mechanism tailored to simulations on Cartesian grids and facilitates
various communication patterns. Individual work steps of a simulation algorithm are specified as sweeps that
are executed on a block-parallel level. The sweep concept defines a structure in which callable objects (i. e.
kernels) implemented in the modules can be specified at compile time. By means of dynamic application
switches, specific kernels tailored to a given computer architecture or implementing a desired model variant,
can be selected at run time. For time-dependent simulations, the sweeps are organized in a timeloop that
specifies the order in which the individual work steps are executed at each time step. To facilitate iterative
solvers, sweeps can also be nested to repeatedly perform a grid traversal until a termination criterion is met.
The boundary condition concept for handling multiple physical fields, numerical methods, and governing
equations, introduced in [31], is applied in this article for moving obstacles with electric BCs. This concept
relies on flags to indicate for each boundary cell the kind of boundary treatment that to is be performed,
18

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

19

with an individual BC flag for each condition. Moreover, cells adjacent to a boundary are indicated with a
nearBC flag and non-boundary cells with a nonBC flag. Individual sets of these flags are defined for each
governing equation. Due to specific LBM requirements, the boundary handling is performed such that the
BCs are fulfilled when a boundary cell is accessed from a nearBC cell in the susequent sweep. The abstract
boundary handling is implemented in the bc module and provides functionality for all BCs are handled
either as direct or direction-dependent treatment. In direction-dependent treatment the BC value is set
at a boundary cell dependent on the value at a neighboring cell, whereas in direct BC treatment the BC
value is directly set at a boundary cell [31]. The actual boundary handling functionality is implemented in
corresponding BC classes whose functions are executed when the associated BC flag is found.
The parameterization concept for multiphysics simulations introduced in [76] is based on the conversion
of physical parameters to lattice units before the simulation, as described in Sec. 4.4. This approach ensures
consistent parameters in all waLBerla modules, independent of the underlying physics. Individual modules
can therefore be developed independently w.r.t. the common lattice unit system.
Simulation parameters and BCs are typically provided to waLBerla through an input file. To ensure
a consistent parameter set and a correct mapping of the physical quantities to lattice units, the class
PhysicalCheck has been introduced in waLBerla in [111]. This class checks the simulation parameter set
specified in the input file for completeness and physical validity and converts the parameters to lattice units.
Since the quantities are converted based on the SI system, the unit Ampere is re-defined for PhysicalCheck
according to Eqn. (53) to support simulations including electric effects.
5.2. Overview of waLBerla Modules for Electrophoresis Simulations
In the following an overview of the modules relevant for electrophoresis simulations is given. For fluid
simulations with the LBM, the lbm module implements various kernels for the stream-collide step with the
different collision operators and forcing terms described in Sec. 4.1. The classes for treating the corresponding
BCs are provided by the associated lbm bc module. In the lbm module block-fields of cells are provided for the
PDFs, the velocity, the density, and an external force. The external force field is used for coupling the LBM
to other methods, e. g. via the forces exerted by electric fields on the EDL and the fluid in electrophoresis
simulations. The PDF and velocity field are accessed by moving obstacle module functions for the simulation
of moving particles.
The moving obstacle module facilitates simulations of fluid-particle interactions with the momentum
exchange method by implementing kernels for moving obstacle sweeps and providing the corresponding
data structures. This module furthermore provides setup functions for initializing and connecting the pe
to waLBerla. For the moving obstacle handling, an obstacle-cell relation field is provided that stores for
each lattice cell overlapped by a pe object a pointer to this object. Moreover, from the lbm module the PDF
source field is accessed for the hydrodynamic force computation and the reconstruction of PDFs. In the lbm
velocity field, body velocities are stored and accessed in the moving boundary treatment of the LBM.
For the computation of the electric potential distribution, the lse solver module described in Sec. 5.3 is
employed. This module has been implemented for solving large sparse linear systems of equations as they
arise from the discretization of the electric potential equations (see Sec. 4.3). The corresponding BC classes
are implemented in the pot bc module described in Sec. 5.4. The data structures defined in the lse solver
module, accessed by the application and other modules, include the stencil field representing the system
matrix as well as the scalar fields for the solution and for the RHS.
The electrokin flow module was designed for facilitating electrokinetic flow simulations. This module
provides kernels for coupling the involved methods as well as the setup and parameterization of simulations
such as electrophoresis. The setup includes initializing the stencils and RHS from the lse solver module
according to the finite volume discretization of the Debye-Hückel equation presented in Sec. 4.3. The kernels
for imposing the electric potential BCs on the moving particles and for computing the electrostatic forces
on fluid and particles are described in Sec. 5.5 and Sec. 5.6, respectively. Finally, the coupled algorithm for
electrophoresis provided by the electrokin flow module is presented in Sec. 5.7.

19

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

20

5.3. Solver for Linear Elliptic PDEs with Moving Boundaries
The solver module lse solver in waLBerla for large linear systems of equations has been designed as
an efficient and robust black-box solver on block-structured grids [112, 31, 76]. To specify the problem to
be solved, the application sets up the system matrix, right-hand side, and BC handling.
In the lse solver module, solver sweeps are pre-defined that perform the iterations at the position where
they are added to the timeloop. For a given simulation setup, the employed solver is selected via the input
file where also the solver parameters and BCs are specified. An iterative solver requires a nested sweep
that is executed until a specific convergence criterion is satisfied. For all implemented solvers, these sweeps
share a common structure. This structure is displayed for the SOR sweep solveTimeVaryingBCSOR for moving
boundaries and linear PDEs such as the Debye-Hückel approximation in Fig. 3 as an activity diagram. For

Adapt stencils
and RHS to BCs

Maximum iterations reached?

yes

no
Compute
residual and norm

Termination criterion met?

yes

no
Update ‘red’
unknowns

Update ‘black’
unknowns
Figure 3: Activity diagram for SOR sweep solveTimeVaryingBCSOR with time varying boundary conditions (BCs),
e. g., due to moving particles

the employed discretization of the electric potential equations, the system matrix is represented by D3Q7
stencils. These stencils comprise an entry for each, the center and the six cardinal directions.
In the setup function for this SOR sweep, the communication for the ghost layer exchange of the solution
field in the initialization phase is set up first. Then the SOR solver sweep is added to the timeloop, and the
kernels for relaxation, communication, and BC treatment are specified as solver sub-sweeps. For parallel
execution, the SOR algorithm is implemented in red-black order.
The filled circle at the top of the diagram in Fig. 3 indicates the starting point of the sweep in the
timeloop. At the beginning of the sweep solveTimeVaryingBCSOR for moving boundaries, the stencils are
constructed to incorporate the BCs according to the present boundary locations. Furthermore, the RHS
is adapted to these BCs before the solver iterations start. For this purpose, a sub-sweep with a kernel
for re-setting the stencils and RHS is executed before the iteration loop, followed by the BC treatment
functions adaptStencilsBC and adaptRHSBC described in Sec. 5.4. Then the standard parallel Red-Black SOR
sweep is performed until the termination criterion is met. This sweep comprises a sub-sweep for computing
the residual and its L2 -norm for the termination criterion, as well as two solver sub-sweeps for the SOR
update of the ‘red’ and the ‘black’ unknowns, respectively. In these sub-sweeps, the quasi-constant stencil
optimization technique introduced in [31] is employed. Based on the residual L2 -norm, the termination
20

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

21

criterion for the simulations performed in this article is provided as residual reduction factor RESRF w.r.t.
the initial norm of the simulation.
For multigrid solvers as applied in [31] to charged particle simulations in absence of ions in the fluid, the
red-black update sub-sweeps in Fig. 3 are replaced by solver sub-sweeps of a V-cycle. To apply the SOR
sweep to varying stencils of linearized PDEs, such as the symmetric Poisson-Boltzmann equation, only the
sub-sweep for adapting the stencils and RHS to the BCs is performed at the beginning of each iteration
instead of once before the iterations start (see [76]).
5.4. Boundary Condition Handling for Solver Module
The implicit BC handling used and initiated by the solver module has been introduced in [31]. This
boundary handling is based on incorporating the BCs into the stencils and right-hand side of the finite
volume discretization. That way, at an iterative update of a near-boundary value, the method implicitly
uses the new values for the BCs. For Dirichlet boundaries, the boundary values are linearly extrapolated to
the boundary cell and for Neumann BCs the boundary values are approximated by central differences. For
both the stencils and the right-hand side, a direction-dependent BC treatment is used.
The functions for this BC treatment are implemented in the pot bc module. This module employs own
nonBC and nearBC flags for the BC handling of scalar potentials. Moreover, for each BC class in this
module an associated BC flag is defined. For the employed cell-centered discretization, the module contains
one class for each, Neumann and Dirichlet BCs.
For incorporating the BCs into the stencils the kernel adaptStencilsBC is implemented. This kernel
iterates over all lattice cells to find scalar potential nearBC cells. At each cell with nearBC flag, the kernel
employs the D3Q7 stencil directions to iterate over the neighboring cells. In directions of a cell with scalar
potential BC flag, the stencil entry of the nearBC cell, associated with the direction of the BC flag, is
adapted accordingly.
The function adaptRHSBC employs the standard boundary handling of the bc module to invoke the pot bc
kernels for adapting the RHS to the BCs. To this end, the direction-dependent BC treatment kernels in
the corresponding BC classes implement the RHS adaption depending on the BC value. The BC value
is specified in the input file for static BCs, or in a previous time step for BCs of moving particles (see
Sec. 5.5). To facilitate such complex boundaries, the BC classes store the BC values and the corresponding
boundary cell ranges. The latter are stored in a memory-efficient way either as cell intervals or as cell vectors.
Moreover, to allow the computation of scalar potential gradients directly from the solution field, the BC
values are set in this field at boundary cells when the RHS is adapted. From these BC values and from the
solution at the nearBC cell, the value at the boundary cell required for the gradient can be extrapolated.
5.5. Electric Potential Boundary Condition Handling for Moving Particles
Prior to the EDL potential distribution computation by the lse solver module, uniform ζ-potentials
or surface charge densities are imposed at the moving particles by means of scalar potential BCs. These
electrical surface properties are specified in the input file for different particle types with a common uid
(unique identifier) defined in the pe. To this end, the sweep function setPotBC_ChargParticles is implemented
in the electrokin flow module that maps the charged particles onto the lattice and sets the electric potential
BC values and the associated BC handling flags at the corresponding cells. The function first overwrites the
values of all cells in the RHS field with zero to remove the values from the previous BC treatment. Then
the mapping is performed for all movable rigid bodies located in the subdomain of each process. For each
particle the mapping is conducted in an extended axis-aligned bounding box that surrounds this rigid body
and the associated nearBC flags. The mapping is realized in three steps:
1) First, all scalar potential nearBC and BC flags from the previous time step are removed, and the
scalar potential nonBC flags are set. Moreover, the previous BC values and the associated cells in the
BC class instances are removed.
2) Then, for particles with prescribed electric BCs, the BC handling for the lse solver module (see Sec. 5.4)
is prepared. For each rigid body with a uid for which a surface property is specified in the input file,
the associated BC is obtained. The cells overlapped by this particle are gathered and are added
21

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

22

together with the BC value to the corresponding BC class instance. Moreover, the BC flag is set at
the overlapped cells.
3) Finally, the nearBC flag is set at all cells adjacent to a BC cell.
Each step is performed for all bodies on a process before the next step begins, to prevent that for particles
with overlapping bounding boxes the flags from a previous step are overwritten.
5.6. Computing Electrostatic Forces on Fluid and Particles
The electric forces acting on the ions in the fluid are incorporated into the incompressible NavierStokes equation by the body force term f~b = −ρe (ψ) ∇ (ϕ + ψ), which coincides with the corresponding
term employed in Henry’s solution (see Sec. 3). Due to the linear superposition of the electric potential
~ ext = −∇ϕ is
components, the gradient is applied to both components separately. Since the applied field E
given, only the EDL potential gradient must be computed.
For the computation of the electric field due to the EDL, the electrokin flow module provides a kernel
that performs the gradient computation at all scalar potential nonBC cells. The gradient of the electric
potential is computed, as previously introduced in [31], by means of finite differences that provide O(dx2 )
accuracy. Where possible, an isotropy-preserving D3Q19 stencil is used (cf. [113]) instead of a D3Q7 stencil.
With the LBM D3Q19 stencil, the gradient can be computed using wq -weighted differences of neighboring
values in 18 directions ~eq as
19
~eq
1 X
wq ψ(~xb + ~eq ) · 2 .
(54)
∇ψ(~xb ) ≈
w1 q=2
dx
At nearBC cells the D3Q7 stencil is applied to compute the gradient of ψ from the BC values stored at
particle and boundary cells in the BC treatment (see Sec. 5.4). The obtained electric field is stored in a field
of cells that is accessed in the body force computation.
For the computation of the electric body force and of the electrostatic force exerted on the particles, a
further kernel is implemented in the electrokin flow module:
The kernel first iterates on each parallel process over all lattice cells to compute the body force at scalar
potential nonBC cells. This force is computed as product of charge density ρe (ψ) and total electric field
~ total = E
~ ext − ∇ψ. The relation of charge density and EDL potential follows Eqn. (20). The obtained
E
electric body force is written to the external force field of the lbm module that is accessed by the LBM
kernels with forcing. Then the kernel iterates over all non-fixed particles residing on the current parallel
process to compute the electrostatic force. For each of these particles the force is computed from the particle
~ ext and is then added to the particle.
charge and the applied field as F~C = qe E
5.7. Algorithm for Electrophoresis
The overall parallel algorithm for electrophoresis simulations with waLBerla is shown in Alg. 1. After
the setup and initialization phase, the electric BCs for the EDL potential are set at the moving charged
particles by means of setPotBC_ChargParticles at each time step. Then the Debye-Hückel approximation is
solved by means of the SOR sweep solveTimeVaryingBCSOR. The iterations are performed until the specified
termination criterion for the residual is met.
From the obtained EDL potential distribution and the applied field the electric body force exerted on the
fluid is computed, as described in Sec. 5.6. First the kernel computing the electric field caused by the EDL is
applied. Then the kernel for the electric body force computation from the charge density distribution in the
fluid and the total electric field is invoked. This kernel additionally applies the electrostatic force exerted
by the applied field to the particles.
Then the rigid body mapping sweep described in Sec. 4.2 is performed, imposing the particle velocities
for the subsequent LBM sweep. In that parallel sweep, an LBM kernel with fused stream-collide step and
forcing term (see Sec. 4.1) is employed to compute the fluid motion influenced by the moving particles and
by the electrostatic force exerted on ions in the EDL.
After the LBM sweep, the hydrodynamic forces on the particles are computed by the momentum exchange
method. The obtained hydrodynamic force contributions and the electrostatic forces are then aggregated by
22

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

23

Algorithm 1: Electrophoresis Algorithm
foreach time step, do
// solve Debye-Hückel approximation (DHA):
set electric BCs of particles
while residual ≥ tol do
apply SOR iteration to DHA
// couple potential solver and LBM:
begin
compute electric field due to EDL
compute charge density in fluid
compute electric body force
apply electrostatic force to particles
// solve lattice Boltzmann equation with forcing, considering particle velocities:
set velocity BCs of particles
begin
perform fused stream-collide step
// couple potential solver and LBM to pe:
begin
apply hydrodynamic force to particles
pe moves particles depending on forces

the pe. From the resulting forces and torques, the new particle velocities and positions are computed in the
subsequent pe simulation step by the PFFD algorithm [114] that additionally resolves rigid body collisions.
6. Electrophoresis Simulations
In the following, the electric potential computation is validated for a sphere with uniform surface charge
surrounded by an EDL. This sphere is placed in a micro-channel subject to an applied electric field. Moreover
the flow field caused by the electrophoretic motion of the sphere in the micro-channel is visualized, together
with the electric potential and the surrounding ions to qualitatively show the correctness of the simulations.
Then the electrophoretic velocity and the retardation by the counter-ions in the EDL is validated w.r.t.
Henry’s solution for different sphere sizes and values of κR.
The simulation setup for the validation experiments is depicted in Fig. 4. A sphere is placed on the
2
L x/

BC

Fy

BC

C
B

z
x

y

BC
BC

Lz/2

Lz/2

2
L x/

C
B

Ly

Figure 4: Setup for electrophoresis of spheres in square duct with different BCs.

longitudinal axis of a cuboid domain of size Lx × Ly × Lz at an initial position of yinit , and an electrostatic
force in y-direction acts on the sphere. The sphere is suspended in a symmetric aqueous electrolyte solution
23

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

24

kg
of kinematic viscosity νf = 1.00 · 10−6 ms , density ρf = 1.00 · 103 m
3 , permittivity εe = 78.54 · ε0 , and ions
of valence z = 1 at a temperature of T = 293 K.
In all simulations, the EDL thickness λD is in the order of the particle diameter. The electric potential
around the sphere due to the EDL is represented by the Debye-Hückel approximation. For this approximation, analytical solutions are known for the potential distribution and for the electrophoretic velocity.
The LBM is employed with TRT operator and second-order forcing term by Luo and re-defined momentum
density (see Sec. 4.1) in all simulations. The linear system of equations resulting from the finite volume
discretization of the Debye-Hückel equation is solved by the SOR method that is sufficient for the quickly
decaying electric potential due to the counter-ions in the EDL. For the SOR, a relaxation parameter of
ωSOR = 1.7 is applied in all simulations.
2

6.1. Electric Potential in an EDL around a Sphere
To validate the computation of the EDL potential ψ around a charged particle, a sphere of radius
R with uniform ζ-potential is simulated in a large domain. The analytical solution of the Debye-Hückel
equation (21) representing the EDL potential around the sphere is given by Eqn. (23). For the validation,
a spherical particle of radius RL = 12 with initial position yinit = 64 dx is chosen and a domain size of
128 dx × 256 dx × 128 dx. The lattice spacing dx is displayed in Tab. 3 with the simulation parameters as
employed related to the electric potential. The ζ-potential is chosen sufficiently low to approximate the
Table 3: Parameters for sphere with uniform surface charge in micro-channel.

c∞ / mol
l

ζ/mV
−10.0

5.00 · 10

1
κ/ m

−6

7.41 · 10

dx/m
6

10 · 10

λD,L

−9

13.49

Poisson-Boltzmann equation by the Debye-Hückel approximation. To obtain the displayed values of κ and
of the characteristic EDL thickness λD , a symmetric aqueous electrolyte solution with ions of the valence
z and the bulk concentration c∞ shown in Tab. 3 is simulated. The EDL thickness is greater than the 12
lattice sites required for a sufficient resolution as observed in the electro-osmotic flow simulations in [76].
Despite its quick decay, the analytical solution of the electric potential at the domain boundaries differs from zero. Thus, the values of ψ at these boundaries are set to the analytical solution by means of
Dirichlet conditions. To solve the Debye-Hückel equation subject to these BCs, a residual reduction factor
of RESRF = 2 · 10−7 is employed as termination criterion for the SOR method.
The analytical (ψ) and numerical (ψ ∗ ) solution at the initial particle position are depicted in Fig. 5 along

ψ/V

0

·10−3

Analytical
solution

−5

−10

Numerical
solution

0

20

40

60
xL

80

100

120

Figure 5: Analytical and numerical solution for EDL potential of sphere with uniform surface charge.

24

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

25

a line in x-direction through the sphere center. Both graphs agree very well, showing the correctness of the
finite volume discretization and the applied SOR solver, as well as the boundary handling at the particle
surface. Inside the insulating particle, the electric potential is not computed explicitly, due to the uniform
surface potential and the resulting symmetric distribution of ψ in the sphere. The electrostatic force needed
to compute the particle motion is computed directly from the applied field and the particle charge, instead
of the electric potential gradient as in [31].
6.2. Electrophoresis of a Sphere in a Micro-Channel
The application of the external electric field in the micro-channel setup described in Sec. 6.1 gives rise to
an electrophoretic motion of the sphere. For the simulation of this scenario, the parameters listed in Tab. 4
are employed in addition to the parameters in Tab. 3. The particle is suspended in an electrolyte solution
Table 4: Parameters for electrophoretic motion of sphere in micro-channel.
V
Ey / m

qs /A s
−18

−19.9 · 10

−4.7 · 10

dt/s
6

−12

200 · 10

UEP,L
4.49 · 10−3

with the parameters introduced at the beginning of Sec. 6. From the ζ-potential and the parameters in
Tab. 3, the sphere’s charge qs in Tab. 4 is obtained from the surface charge density σs according to the ζ–σs
relationship (32), multiplied by the surface area of the sphere, as qs = 4πR2 σs [81].
The high electric field strength Ey in Tab. 4 is chosen to keep the number of simulation time steps at a
minimum. Due to the applied field and the resulting electrostatic force FC,y = 933 · 10−12 N the particle
moves in y-direction, retarded by the channel walls and the opposing force on the EDL. For the chosen
parameters, the terminal particle speed of UEP = 224.5 mm
s is obtained for free space according to Henry’s
solution (see Eqn. (29)), corresponding to a particle Reynolds number of Rep,d = 0.054.
In the simulation, periodic BCs are applied in y-direction. At all other walls, no-slip conditions are
applied for the LBM and homogeneous Neumann conditions for the electric potential. At each time step,
the Debye-Hückel equation is solved by SOR with the termination criterion from Sec. 6.1. The LBM is
run with τ = 6.5, resulting in the time increment dt listed in Tab. 4 for the chosen viscosity νf and dx
(cf. Eqn. (36)). With dt and dx, the electrophoretic speed in lattice units UEP,L attains the values given
in Tab. 4. Gravitational effects are neglected in the simulation to ensure that the particle motion is driven
kg
solely by electric forces. Thus, the employed density ρp = 1195 m
3 of the particle only has an impact on the
time required to reach steady-state.
In Fig. 6, the results of the electrophoresis simulation are visualized at different time steps. The EDL
potential ψ in the y-z plane through the domain center is displayed, together with semi-transparent equipotential surfaces representing the excess counter-ions in the EDL. The flow field around the moving sphere is
visualized in the x-z plane through the domain center. Arrows of uniform length indicate the flow direction,
while the velocity magnitude is represented by the shown color-scale and by twelve white isosurface contour
lines with logarithmic intervals in the range of 13.67 · 10−3 ms to 6.08 · 10−6 ms .
The flow field around the moving particle shown in Fig. 6(a) is fully developed after 5001 time steps, and
the particle has already attained its terminal velocity. Due to the periodicity of the domain, a channel flow in
axial direction has emerged from the particle motion, as indicated by the contour lines. Moreover, a vortex
has formed between the sphere and the surrounding no-slip walls. The flow field moves with the particle
that translates along the channel centerline (see Fig. 6(b)). Because of the equilibrium representation of the
EDL, the distribution of ψ is at all time steps symmetric w.r.t. the sphere center, almost up to the boundary.
6.3. Validation of Electrophoretic Motion of a Sphere
To quantitatively validate the overall electrophoresis simulation algorithm, the electrophoretic velocity
of spheres with uniform surface charge is compared to Henry’s solution Eqn. (29) for a spherical particle in
an unbounded electrolyte solution.
25

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

26

(a) Results after 5001 time steps.

(b) Results after 30 001 time steps.
Figure 6: Electrophoresis of spherical particle in micro-channel with insulating no-slip walls. Visualization of flow field in x-y
plane, EDL potential in y-z plane, and ion charge distribution around charged particle.

In the simulation experiments, a sphere with uniform ζ-potential is moving under the influence of an
~ ext in a large domain filled with an electrolyte solution. The simulations are
applied field of strength E
performed for different values of κR by varying the sphere radius while keeping the EDL thickness constant.
For validation, the relative deviation ∆r U of the obtained terminal sphere velocity from the theoretical
velocity in Eqn. (29) is evaluated.
The simulation parameters are chosen such that the electrophoretic motion is in the Stokes regime, and
thus the nonlinear inertial term of Navier-Stokes equation can be neglected. Moreover, the EDL potential
distribution is governed by the linear Debye-Hückel approximation. Therefore, the superposition principle is
assumed to hold. Thus, from the obtained relative deviation from the analytical solution in an unbounded
domain ∆r U , the relative deviation due to wall effects and volume mapping errors ∆r UStokes will be subtracted. These relative deviations were examined in [76] for several domain sizes and sphere radii. For no-slip
BCs, the wall effect was shown to comply with analytical and experimental results. In the experiments,
domain sizes are used for which the relative deviations ∆r UStokes for free-slip BCs are close to −3 %.
For all simulations, the parameters in Tab. 5 are used. In each experiment, a sphere is suspended in a
symmetric aqueous electrolyte solution with the parameters introduced at the beginning of Sec. 6 and a bulk
concentration of c∞ = 1.60 · 10−5 mol
l . These parameters result in the Debye-Hückel parameter κ shown in
Tab. 5 and, together with the displayed lattice spacing dx, in the EDL thickness λD of approximately 15
26

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

27

Table 5: Parameters for electrophoretic velocity validation w.r.t. Henry’s solution.

ζ/mV

1
κ/ m

V
Ey / m

dx/m

λD,L

10.0

13.3 · 106

99.0 · 106

5.00 · 10−9

15.08

lattice sites. The ζ-potential has the same absolute value as in the electrophoresis simulations in a microchannel from Sec. 6.2. Moreover, as in Sec. 6.2, a high electric field strength Ey is chosen to keep the number
of simulation time steps low. The resulting electrostatic forces FC = qs Ey for the different sphere radii
are displayed in Tab. 6, together with the electrophoretic charges qs on the spheres. These surface charges
associated with the ζ-potential are obtained from the general ζ–σs relationship (32) as in Sec. 6.2. For all
simulation parameters, the maximum relative deviation of σs for the general relationship from the exact
value for the Debye-Hückel approximation is below 0.2 %. These charges are applied as particle charges to
drive the spheres by the electrostatic force due to the electric field.
For the different sphere radii and the associated values of κR, the electrophoretic velocities according
to Henry’s solution are displayed in Tab. 6. These velocities correspond to particle Reynolds numbers Rep,d
from 0.018 to 0.057 for the particle diameters of 40 nm to 120 nm.
EM
is introduced.
To quantify the retardation by the opposing force on the EDL, the variable EPRet = UEPU−U
EM
This variable represents the relative deviation of the electrophoretic velocity of a particle with charge qs in
presence of the electric double layer from the migration velocity UEM of a particle with the same charge in
absence of surrounding ions. For the examined sphere radii, this retardation is in the range of 20 % to 42 %.
Table 6: Electrophoresis parameters and domain sizes dependent on sphere size. Shown are the sphere radii RL in lattice units,
the corresponding charges qs , and the electrostatic forces FC . For relation of sphere radius to EDL thickness κR, theoretical
electrophoretic velocities UEP , Reynolds numbers Rep,d , and electrophoretic retardation parameter EPRet are given. Listed in
lower part are domain sizes per dimension Lx,y,z , initial sphere positions yinit , process numbers per dimension #procx,y,z , and
relative deviations of sphere velocities in free-slip domain from Stokes velocity ∆r UStokes .

4

6

8

9

12

qs /(10
A s)
FC /(10−12 N)

2.21
219

3.67
363

5.36
530

6.29
622

9.43
934

κR
UEP /(10−3
Rep,d
EPRet /%

0.265
461
0.018
-20.6

0.398
464
0.028
-27.7

0.530
466
0.037
-33.6

0.597
468
0.042
-36.2

0.796
471
0.057
-42.8

864
1248
392

1280
1536
588

1632
2048
784

1632
2048
882

1632
2048
1176

8 × 16 × 8

8 × 16 × 16

16 × 16 × 32

16 × 16 × 32

16 × 16 × 32

-3.04

-2.65

-2.98

-2.89

-3.00

RL
−18

m
s )

Lx,z /dx
Ly /dx
yinit /dx
#procx,y,z
∆r UStokes /%

Then the domain sizes for which the relative deviation from Stokes velocity is about ∆r UStokes = −3 %
are listed in Tab. 6, together with the initial position yinit in movement direction that correspond to 98 × RL .
The LBM with TRT collision operator is employed with the relaxation time τ = 6, resulting in the time
increment dt = 45.8 · 10−12 s. As in the micro-channel electrophoresis simulations, gravitational are effects
kg
neglected, and the insulating particles have a density of ρp = 1195 m
3.
For the electric potential, homogeneous Neumann BCs are applied at the walls in x- and z-direction as
27

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

28

= 0. To improve convergence [31], homogeneous Dirichlet BCs are applied at the walls
in y-direction as Φ Γ ∪Γ = 0 V. The particle is located at sufficient distance from these walls for the
N
S
EDL potential to decay to approximately zero, and therefore these BCs do not affect the electric potential
distribution. As termination criterion of the SOR solver for the Debye-Hückel equation, a residual reduction
factor of RESRF = 1 · 10−6 is used.
The parallel simulations are run on SuperMUC1 of the Leibniz Supercomputing Centre LRZ in Garching
(Germany) on 64 to 512 nodes with the numbers of processes listed in Tab. 6. Within the execution time
of 37 h to 48 h, a number of 70 296 to 74 000 time steps are performed, and the spheres cover a distance
of about 300 dx. The high numbers of time steps are chosen to ensure that the spheres reach steady-state
motion, and that large numbers of sampling values are available to compute the terminal velocities.
The sphere velocities are sampled every 20 time steps during the simulation. In the second half of the
simulation, the terminal particle velocity is reached. Thus, the mean particle velocity U ∗ and the velocity
U ∗ −U ∗
fluctuations δU = maxU ∗ min due to volume mapping effects are computed from the last 50 % output values.
In the considered range of time steps additionally the number of time steps between two SOR calls and the
number of SOR iterations is monitored. The average number of time steps between two SOR calls decreases
from ∆TSSOR = 24 to ∆TSSOR = 3 as RL increases from 4 to 12. Likewise, the average number of SOR
iterations per solver call decreases from 451 iterations to 198 iterations for the respective sphere radii.
The obtained terminal velocity in lattice units UL∗ and the fluctuations are displayed in Tab. 7. As expected,
∂Φ
∂~
n ΓW ∪ΓE ∪ΓB ∪ΓT

Table 7: Simulation results of electrophoretic velocity validation for different sphere sizes. Shown are the theoretical velocities
∗ and fluctuations δ , the relative deviations ∆ U of
UEP,L for unbounded domains in lattice units, the obtained velocities UL
r
U
∗ from U
UL
EP,L , and the relative deviations ∆r UEP corrected by hydrodynamic wall and mapping effects.

4

6

8

9

12

UEP,L /10

4.227

4.249

4.274

4.286

4.320

UL∗ /10−3

4.119
2.52
-2.56
0.5

4.144
1.52
-2.48
0.2

4.120
0.83
-3.59
-0.6

4.135
0.68
-3.51
-0.6

4.151
0.29
-3.91
-0.9

RL
−3

δU /%
∆r U/%
∆r UEP /%

the fluctuations decrease with increasing sphere resolution. Moreover, the relative deviation of the obtained
velocity U ∗ from the theoretical electrophoretic velocity UEP given by ∆r U = (U ∗ − UEP )/UEP is listed
in Tab. 7. For all examined sphere radii, the obtained velocities in the confined domain are by 2.5 % to
3.9 % lower than the theoretical values of UEP for a particle in an unbounded electrolyte solution. The
effect of the confinement on the particle velocity is deducted by subtracting the relative deviation from
Stokes velocity ∆r UStokes in Tab. 6 for the corresponding domain sizes from the relative electrophoretic
velocity deviation ∆r U obtained in the electrophoresis simulation. From the resulting relative deviations
∆r UEP = ∆r U − ∆r UStokes , the inaccuracies due to electric effects in the electrophoresis simulation are
assessed. As can be seen from the values of ∆r UEP in Tab. 7, the simulations results agree with the theoretical
values with relative deviations of less than 1 %.
7. Conclusion
In this article, a coupled multiphysics algorithm for parallel simulations of the electrophoretic motion
of geometrically resolved particles in electrolyte solutions is presented. The physical effects are simulated
by means of the lattice Boltzmann method for fluid dynamics, a physics engine for rigid body dynamics,
and a scalar iterative solver for electric potentials. These components are integrated into the parallel
software framework waLBerla to simulate the electrophoretic motion of fully resolved charged particles in
1 www.lrz.de/services/compute/supermuc/

28

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

29

microfluidic flow. The simulations include fluid-particle and electrostatic particle interactions. Additionally
electric effects on ions around the homogeneously charged particles are recovered.
The current work is an extension of [31], where the electrical migration of charged particles without
ions in the fluid was validated and excellent parallel performance was shown for more than seven millions
of interacting charged particles. In the present article, the opposite net charge in the electric double layer
(EDL) around the charged particles due to ions in the fluid is considered, together with its effect on the fluid
motion that counteracts particle motion. To this end, the electric potential distribution in the fluid due to
the EDLs is computed that causes an electric body force on the fluid. This quasi-equilibrium distribution
recovers the motion of ions in the fluid along with the charged particles while neglecting EDL distortion.
The overall electrophoresis algorithm is introduced and an overview of the coupled functionality implemented in the involved waLBerla modules is given. For the simulations, a solver sweep for time-varying
boundary conditions has been developed that is presented here for the parallel SOR method employed to
solve the EDL potential equation. Based on the multiphysics boundary handling concept [31] an efficient
parallel algorithm is implemented to impose electric potential boundary conditions on the moving particles. These methods can also be employed for other governing equations with spatially varying boundary
conditions that model physical effects different from electric fields. The presented parallel electrophoresis simulations are also facilitated by a joint parameterization concept for the different coupled governing
equations and numerical methods implemented in waLBerla. This concept is based on lattice Boltzmann
requirements and is applicable and extensible to further multiphysics simulations.
For the electrophoresis simulations in this article, the electric potential in the double layer is shown to
coincide with analytical solutions. The obtained terminal electrophoretic velocities comply with analytical
solutions for different proportions of the particle radii to double layer thickness. These validation results
verify the correctness of the implementation and the coupling of the different methods. Moreover, the
observed relative errors in the modeling of electric effects are below 1 %. The retardation effect caused
by the presence of the EDL is shown to be significant for the examined sphere radii, reducing the sphere
velocity up to 42 %. For the electrophoretic motion in a micro-channel, the flow field and the electric potential
distribution are visualized, including the ion charge distribution in the EDL surrounding the particle.
The presented algorithm can be applied to find design parameters in industrial and medical applications,
e. g., for optimal separation efficiency of charged biological particles in lab-on-a-chip devices, depending on
fluid, particle, and electrolyte properties. Our algorithms were shown to correctly recover fluid-particle
interactions for elongated particles in [32]. These future simulations may therefore include suspended,
possibly charged particles of various shapes including spheres, spherocylinders, and particles of more complex
shapes, e. g., to represent different biological particles. Also pairwise van-der-Waals forces can be added
easily, to facilitate simulations of electrophoretic deposition in material science applications.
The electrophoresis algorithm introduced here is well suited for massively parallel simulations. In the
current implementation of this algorithm, the EDL thickness is restricted to values in the order of the particle
radius. Therefore, adaptive lattice refinement as in [23] may be employed to allow for thinner double layers
relative to the particle size. For the incorporation of transient effects in the simulations including EDLs,
the link-flux method implemented into waLBerla in [115] and [41] may be employed. This method was
extended in [41] to simulate electrophoresis, enabling the simulation of non-equilibrium ion distributions in
the EDL. Due to the higher computational complexity of the link-flux method compared to the equilibrium
approach in this article, the maximum number of particles will be lower than in our approach.
Acknowledgements
The authors are grateful to the LRZ for providing the computational resources on SuperMUC.
References
[1] K. J. Ptasinski, P. J. A. M. Kerkhof, Electric Field Driven Separations: Phenomena and Applications, Sep. Sci. Technol.
27 (8-9) (1992) 995–1021. doi:10.1080/01496399208019021.

29

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

30

[2] S. Zhang, R. Tan, K. Neoh, C. Tien, Electrofiltration of Aqueous Suspensions, J. Colloid Interface Sci. 228 (2) (2000)
393–404. doi:10.1006/jcis.2000.6966.
[3] Y.-H. Weng, K.-C. Li, L. H. Chaung-Hsieh, C. Huang, Removal of humic substances (HS) from water by electromicrofiltration (EMF), Water Res. 40 (9) (2006) 1783–1794. doi:10.1016/j.watres.2006.02.028.
[4] A. Mahmoud, J. Olivier, J. Vaxelaire, A. Hoadley, Electrical field: A historical review of its application and contributions
in wastewater sludge dewatering, Water Res. 44 (8) (2010) 2381–2407. doi:10.1016/j.watres.2010.01.033.
[5] L. Besra, M. Liu, A review on fundamentals and applications of electrophoretic deposition (EPD), Prog. Mater. Sci.
52 (1) (2007) 1 – 61. doi:10.1016/j.pmatsci.2006.07.001.
[6] I. Zhitomirsky, Cathodic electrodeposition of ceramic and organoceramic materials. Fundamental aspects, Adv. Colloid
Interface Sci. 97 (1–3) (2002) 279–317. doi:10.1016/S0001-8686(01)00068-9.
[7] P. Sarkar, D. De, T. Uchikochi, L. Besra, Electrophoretic Deposition (EPD): Fundamentals and Novel Applications in
Fabrication of Advanced Ceramic Microstructures, Springer, 2012, Ch. 5, pp. 181–215. doi:10.1007/978-1-4419-9730-2_
5.
[8] Y. Kang, D. Li, Electrokinetic motion of particles and cells in microchannels, Microfluid Nanofluid. 6 (4) (2009) 431–460.
doi:10.1007/s10404-009-0408-7.
[9] A. A. S. Bhagat, H. Bow, H. W. Hou, S. J. Tan, J. Han, C. T. Lim, Microfluidics for cell separation, Med. Biol. Eng.
Comput. 48 (10) (2010) 999–1014. doi:10.1007/s11517-010-0611-4.
[10] D. R. Gossett, W. M. Weaver, A. J. Mach, S. C. Hur, H. T. K. Tse, W. Lee, H. Amini, D. Di Carlo, Label-free
cell separation and sorting in microfluidic systems, Anal. Bioanal. Chem. 397 (8) (2010) 3249–3267. doi:10.1007/
s00216-010-3721-9.
[11] N. Pamme, Continuous flow separations in microfluidic devices, Lab Chip 7 (2007) 1644–1659. doi:10.1039/B712784G.
[12] D. G. Hert, C. P. Fredlake, A. E. Barron, Advantages and limitations of next-generation sequencing technologies: A
comparison of electrophoresis and non-electrophoresis methods, Electrophoresis 29 (23) (2008) 4618–4626. doi:10.1002/
elps.200800456.
[13] G. W. Slater, C. Holm, M. V. Chubynsky, H. W. de Haan, A. Dube, K. Grass, O. A. Hickey, C. Kingsburry, D. Sean,
T. N. Shendruk, L. Zhan, Modeling the separation of macromolecules: A review of current computer simulation methods,
Electrophoresis 30 (5) (2009) 792–818. doi:10.1002/elps.200800673.
[14] F. Keller, H. Nirschl, W. Dörfler, E. Woldt, Efficient numerical simulation and optimization in electrophoretic deposition
processes, J. Eur. Ceram. Soc. 35 (9) (2015) 2619–2630. doi:10.1016/j.jeurceramsoc.2015.02.031.
[15] D. Sheehan, Physical Biochemistry: Principles and Applications, John Wiley & Sons, 2013.
[16] R. F. Probstein, Physicochemical Hydrodynamics: An Introduction, 2nd Edition, Butterworths series in chemical engineering, Butterworth Publishers, 1989.
[17] H. C. Chang, L. Y. Yeo, Electrokinetically Driven Microfluidics and Nanofluidics, Cambridge Univ. Press, 2009.
[18] O. Stern, ZUR THEORIE DER ELEKTROLYTISCHEN DOPPELSCHICHT, Z. Elektrochem. 30 (21-22) (1924) 508–
516. doi:10.1002/bbpc.192400182.
[19] T. Preclik, U. Rüde, Ultrascale Simulations of non-smooth granular dynamics, Comp. Part. Mech. 2 (2) (2015) 173–196.
doi:10.1007/s40571-015-0047-6.
[20] K. Iglberger, U. Rüde, Massively Parallel Rigid Body Dynamics Simulation, Comp. Sci. Res. Dev. 23 (3-4) (2009) 159–167.
doi:10.1007/s00450-009-0066-8.
[21] C. Feichtinger, S. Donath, H. Köstler, J. Götz, U. Rüde, WaLBerla: HPC software design for computational engineering
simulations, J. Comput. Sci. 2 (2) (2011) 105–112. doi:10.1016/j.jocs.2011.01.004.
[22] C. Godenschwager, F. Schornbaum, M. Bauer, H. Köstler, U. Rüde, A Framework for Hybrid Parallel Flow Simulations
with a Trillion Cells in Complex Geometries, in: Proc. Int. Conf. on High Performance Computing, Networking, Storage
and Analysis, SC ’13, ACM, 2013, pp. 35:1–35:12. doi:10.1145/2503210.2503273.
[23] F. Schornbaum, U. Rüde, Massively Parallel Algorithms for the Lattice Boltzmann Method on NonUniform Grids, SIAM
J. Sci. Comput. 38 (2) (2016) C96–C126. arXiv:1508.07982, doi:10.1137/15M1035240.
[24] S. Chen, G. D. Doolen, Lattice Boltzmann Method for Fluid Flows, Annu. Rev. Fluid Mech. 30 (1) (1998) 329–364.
doi:10.1146/annurev.fluid.30.1.329.
[25] C. K. Aidun, J. R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (1) (2010) 439–472.
doi:10.1146/annurev-fluid-121108-145519.
[26] A. J. C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical
foundation, J. Fluid Mech. 271 (1994) 285–309. doi:10.1017/S0022112094001771.
[27] N. Q. Nguyen, A. J. C. Ladd, Lubrication corrections for lattice-Boltzmann simulations of particle suspensions, Phys.
Rev. E 66 (4) (2002) 046708. doi:10.1103/PhysRevE.66.046708.
[28] J. Götz, K. Iglberger, M. Stürmer, U. Rüde, Direct Numerical Simulation of Particulate Flows on 294912 Processor
Cores, in: Proc. 2010 ACM/IEEE Proc. Int. Conf. for High Performance Computing, Networking, Storage and Analysis,
SC ’10, IEEE, 2010, pp. 1–11. doi:10.1109/SC.2010.20.
[29] I. Ginzburg, J.-P. Carlier, C. Kao, Lattice Boltzmann approach to Richards’ equation, in: W. G. G. C. T. Miller,
M. W. Farthing, G. F. Pinder (Eds.), Computational Methods in Water Resources, Vol. 55 of Developments in Water
Science, Elsevier, 2004, pp. 583–595. doi:10.1016/S0167-5648(04)80083-2.
[30] I. Ginzburg, F. Verhaeghe, D. d’Humières, Two-relaxation-time lattice Boltzmann scheme: About parametrization,
velocity, pressure and mixed boundary conditions, Commun. Comput. Phys. 3 (2) (2008) 427–478.
[31] D. Bartuschat, U. Rüde, Parallel Multiphysics Simulations of Charged Particles in Microfluidic Flows, J. Comput. Sci.
8 (0) (2015) 1–19. doi:10.1016/j.jocs.2015.02.006.
[32] D. Bartuschat, E. Fischermeier, K. Gustavsson, U. Rüde, Two Computational Models for Simulating the Tumbling

30

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

31

Motion of Elongated Particles in Fluids, Comput. Fluids 127 (2016) 17–35. doi:10.1016/j.compfluid.2015.12.010.
[33] S. Jomeh, M. Hoorfar, Study of the effect of electric field and electroneutrality on transport of biomolecules in microreactors, Microfluid Nanofluid. 12 (1) (2012) 279–294. doi:10.1007/s10404-011-0871-9.
[34] P. Kler, C. Berli, F. Guarnieri, Modeling and high performance simulation of electrophoretic techniques in microfluidic
chips, Microfluid Nanofluid. 10 (1) (2011) 187–198. doi:10.1007/s10404-010-0660-x.
[35] M. Chau, T. Garcia, P. Spiteri, Asynchronous grid computing for the simulation of the 3D electrophoresis coupled
problem, Adv. Eng. Softw. 60–61 (2013) 111–121. doi:10.1016/j.advengsoft.2012.11.010.
[36] M. Chau, P. Spiteri, H. C. Boisson, Parallel numerical simulation for the coupled problem of continuous flow electrophoresis, Int. J. Numer. Meth. Fluids 55 (10) (2007) 945–963. doi:10.1002/fld.1502.
[37] B. Giera, L. A. Zepeda-Ruiz, A. J. Pascall, J. D. Kuntz, C. M. Spadaccini, T. H. Weisgraber, Mesoscale particle-based
model of electrophoresis, J. Electrochem. Soc. 162 (11) (2015) D3030–D3035. doi:10.1149/2.0161511jes.
[38] R. G. M. van der Sman, Simulations of confined suspension flow at multiple length scales, Soft Matter 5 (2009) 4376–4387.
doi:10.1039/B915749M.
[39] J. Smiatek, F. Schmid, Mesoscopic Simulations of Electroosmotic Flow and Electrophoresis in Nanochannels, Comput.
Phys. Commun. 182 (9) (2011) 1941 – 1944, computer Physics Communications Special Edition for Conference on
Computational Physics Trondheim, Norway, June 23-26, 2010. doi:10.1016/j.cpc.2010.11.021.
[40] R. Wang, J.-S. Wang, G.-R. Liu, J. Han, Y.-Z. Chen, Simulation of DNA electrophoresis in systems of large number
of solvent particles by coarse-grained hybrid molecular dynamics approach, J. Comput. Chem. 30 (4) (2009) 505–513.
doi:10.1002/jcc.21081.
[41] M. Kuron, G. Rempfer, F. Schornbaum, M. Bauer, C. Godenschwager, C. Holm, J. de Graaf, Moving charged particles
in lattice Boltzmann-based electrokinetics, J. Chem. Phys. 145 (21) (2016) 214102. doi:10.1063/1.4968596.
[42] J. S. Park, D. Saintillan, Direct Numerical Simulations of Electrophoretic Deposition of Charged Colloidal Suspensions,
in: R. C. A. R. Boccaccini, O. Van der Biest, J. Dickerson (Eds.), Key Engineering Materials, Vol. 507, Trans Tech Publ,
2012, pp. 47–51. doi:10.4028/www.scientific.net/KEM.507.47.
[43] J.-P. Hsu, C.-H. Chou, C.-C. Kuo, S. Tseng, R. Wu, Electrophoresis of an arbitrarily oriented toroid in an unbounded
electrolyte solution, Colloids Surf., B 82 (2) (2011) 505–512. doi:10.1016/j.colsurfb.2010.10.009.
[44] J.-P. Hsu, M.-H. Ku, Boundary effect on electrophoresis: finite cylinder in a cylindrical pore, J. Colloid Interface Sci.
283 (2) (2005) 592–600. doi:10.1016/j.jcis.2004.09.004.
[45] J.-P. Hsu, L.-H. Yeh, Electrophoresis of Two Identical Rigid Spheres in a Charged Cylindrical Pore, J. Phys. Chem. B
111 (10) (2007) 2579–2586. doi:10.1021/jp068407z.
[46] J. H. Masliyah, S. Bhattacharjee, Electrokinetic and Colloid Transport Phenomena, John Wiley & Sons, Inc., 2006.
[47] S. Tseng, C.-H. Huang, J.-P. Hsu, Electrophoresis of two spheres: Influence of double layer and van der Waals interactions,
J. Colloid Interface Sci. 451 (2015) 170–176. doi:10.1016/j.jcis.2015.03.060.
[48] R. Schmitz, V. Starchenko, B. Dünweg, Computer simulation of electrokinetics in colloidal systems, Eur. Phys. J. Spec.
Topics 222 (11) (2013) 2873–2880. doi:10.1140/epjst/e2013-02063-2.
[49] M. Baptista, R. Schmitz, B. Dünweg, Simple and robust solver for the Poisson-Boltzmann equation, Phys. Rev. E 80
(2009) 016705. doi:10.1103/PhysRevE.80.016705.
[50] H. H. Hu, N. A. Patankar, M. Y. Zhu, Direct Numerical Simulations of Fluid-Solid Systems Using the Arbitrary
Lagrangian-Eulerian Technique, J. Comput. Phys. 169 (2) (2001) 427–462. doi:10.1006/jcph.2000.6592.
[51] T. J. R. Hughes, W. K. Liu, T. K. Zimmermann, Lagrangian-Eulerian finite element formulation for incompressible
viscous flows, Compu. Meth. Appl. Mech. and Engin. 29 (3) (1981) 329–349. doi:10.1016/0045-7825(81)90049-9.
[52] C. Ye, D. Li, 3-D transient electrophoretic motion of a spherical particle in a T-shaped rectangular microchannel, J.
Colloid Interface Sci. 272 (2) (2004) 480–488. doi:10.1016/j.jcis.2003.11.014.
[53] H. Tanaka, T. Araki, Simulation Method of Colloidal Suspensions with Hydrodynamic Interactions: Fluid Particle
Dynamics, Phys. Rev. Lett. 85 (2000) 1338–1341. doi:10.1103/PhysRevLett.85.1338.
[54] H. Kodama, K. Takeshita, T. Araki, H. Tanaka, Fluid particle dynamics simulation of charged colloidal suspensions, J.
Phys. Condens. Matter 16 (10) (2004) L115–L123. doi:10.1088/0953-8984/16/10/L01.
[55] T. Araki, H. Tanaka, Physical principle for optimizing electrophoretic separation of charged particles, Europhys. Lett.
82 (1) (2008) 18004. doi:10.1209/0295-5075/82/18004.
[56] Y. Nakayama, R. Yamamoto, Simulation method to resolve hydrodynamic interactions in colloidal dispersions, Phys.
Rev. E 71 (2005) 036707. doi:10.1103/PhysRevE.71.036707.
[57] K. Kim, Y. Nakayama, R. Yamamoto, Direct Numerical Simulations of Electrophoresis of Charged Colloids, Phys. Rev.
Lett. 96 (2006) 208302. doi:10.1103/PhysRevLett.96.208302.
[58] C. Shih, R. Yamamoto, Dynamic electrophoresis of charged colloids in an oscillating electric field, Phys. Rev. E 89 (2014)
062317. doi:10.1103/PhysRevE.89.062317.
[59] X. Luo, A. Beskok, G. E. Karniadakis, Modeling Electrokinetic Flows by the Smoothed Profile Method, J. Comput.
Phys. 229 (10) (2010) 3828 – 3847. doi:10.1016/j.jcp.2010.01.030.
[60] C. S. Peskin, Flow patterns around heart valves: A numerical method, J. Comput. Phys. 10 (2) (1972) 252–271. doi:
10.1016/0021-9991(72)90065-4.
[61] R. Mittal, G. Iaccarino, Immersed boundary methods, Annu. Rev. Fluid Mech. 37 (2005) 239–261. doi:10.1146/annurev.
fluid.37.061903.175743.
[62] M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. Comput.
Phys. 209 (2) (2005) 448–476. doi:10.1016/j.jcp.2005.03.017.
[63] S. Kang, Direct simulations on the electrophoretic motion of multiple charged particles using an immersed boundary
method, Comput. Fluids 73 (2013) 10–23. doi:10.1016/j.compfluid.2012.12.005.

31

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

32

[64] V. Lobaskin, B. Dünweg, C. Holm, Electrophoretic mobility of a charged colloidal particle: a computer simulation study,
J. Phys. Condens. Matter 16 (38) (2004) S4063–S4073. doi:10.1088/0953-8984/16/38/021.
[65] T. Y. Molotilin, V. Lobaskin, O. I. Vinogradova, Electrophoresis of Janus particles: A molecular dynamics simulation
study, J. Chem. Phys. 145 (24) (2016) 244704. doi:10.1063/1.4972522.
[66] A. Chatterji, J. Horbach, The role of effective charges in the electrophoresis of highly charged colloids, J. Phys. Condens.
Matter 22 (49) (2010) 494102. doi:10.1088/0953-8984/22/49/494102.
[67] O. A. Hickey, C. Holm, J. Smiatek, Lattice-Boltzmann simulations of the electrophoretic stretching of polyelectrolytes:
The importance of hydrodynamic interactions, J. Chem. Phys. 140 (16) (2014) 164904. doi:10.1063/1.4872366.
[68] F. Capuani, I. Pagonabarraga, D. Frenkel, Discrete solution of the electrokinetic equations, J. Chem. Phys. 121 (2) (2004)
973–986. doi:10.1063/1.1760739.
[69] G. Giupponi, I. Pagonabarraga, Determination of the zeta potential for highly charged colloidal suspensions, Phil. Trans.
R. Soc. A: Mathematical, Physical and Engineering Sciences 369 (2011) 2546–2554. doi:10.1098/rsta.2011.0024.
[70] G. K. Batchelor, An Introduction to Fluid Dynamics, Cambridge Univ. Press, 1979.
[71] M. A. van der Hoef, M. Ye, M. van Sint Annaland, A. T. Andrews, S. Sundaresan, J. A. M. Kuipers, Multi-Scale
Modeling of Gas-Fluidized Beds, in: G. B. Marin (Ed.), Computational Fluid Dynamics, Vol. 31 of Advances in Chemical
Engineering, Elsevier, 2006, pp. 65–149. doi:10.1016/S0065-2377(06)31002-2.
[72] L. D. Landau, E. M. Lifshitz, Fluid Mechanics, Second Edition: Volume 6 (Course of Theoretical Physics), 2nd Edition,
Vol. 6 of Course of theoretical physics / by L. D. Landau and E. M. Lifshitz, Vol. 6, Butterworth-Heinemann, 1987.
[73] P. J. Dellar, Bulk and shear viscosities in lattice Boltzmann equations, Phys. Rev. E 64 (3) (2001) 031203. doi:
10.1103/PhysRevE.64.031203.
[74] J. Happel, H. Brenner, Low Reynolds number hydrodynamics: with special applications to particulate media, Vol. 1,
Springer, 1983.
[75] F. Durst, Grundlagen Der Strömungsmechanik, Springer, 2006.
[76] D. Bartuschat, Direct Numerical Simulation of Particle-Laden Electrokinetic Flows on High Performance Computers,
Ph.D. thesis, University of Erlangen-Nürnberg (2016).
URL https://opus4.kobv.de/opus4-fau/frontdoor/index/index/docId/7298
[77] J. Anderson, Computational Fluid Dynamics, Computational Fluid Dynamics: The Basics with Applications, McGrawHill Education, 1995.
[78] G. G. Stokes, On the effect of the internal friction of fluids on the motion of pendulums, Vol. 9, Pitt Press, 1851.
[79] J. Happel, E. Bart, The settling of a sphere along the axis of a long square duct at low Reynolds’ number, Appl. Sci.
Res. 29 (1) (1974) 241–258. doi:10.1007/BF00384149.
[80] A. Miyamura, S. Iwasaki, T. Ishii, Experimental wall correction factors of single solid spheres in triangular and square
cylinders, and parallel plates, Int. J. Multiph. Flow 7 (1) (1981) 41–46. doi:10.1016/0301-9322(81)90013-6.
[81] R. J. Hunter, Zeta Potential in Colloid Science: Principles and Applications, Academic Press, London, 1981.
[82] P. Debye, E. Hückel, Zur Theorie der Elektrolyte. I. Gefrierpunktserniedrigung und verwandte Erscheinungen, Phys. Z
24 (9) (1923) 185–206.
[83] D. C. Henry, The cataphoresis of suspended particles. Part I. The equation of cataphoresis, Proc. R. Soc. Lond. A
133 (821) (1931) 106–129. doi:10.1098/rspa.1931.0133.
[84] H. Ohshima, A Simple Expression for Henry’s Function for the Retardation Effect in Electrophoresis of Spherical Colloidal
Particles, J. Colloid Interface Sci. 168 (1) (1994) 269–271. doi:10.1006/jcis.1994.1419.
[85] H. Ohshima (Ed.), Theory of Colloid and Interfacial Electric Phenomena, Vol. 12 of Interface Science and Technology,
Elsevier, 2006. doi:10.1016/S1573-4285(06)80022-0.
[86] D. J. Shaw, Introduction to Colloid and Surface Chemistry, 4th Edition, Butterworth-Heinemann, 1992.
[87] H. Ohshima, T. W. Healy, L. R. White, Accurate analytic expressions for the surface charge density/surface potential
relationship and double-layer potential distribution for a spherical colloidal particle, J. Colloid Interface Sci. 90 (1) (1982)
17–26. doi:10.1016/0021-9797(82)90393-9.
[88] A. L. Loeb, J. T. G. Overbeek, P. H. Wiersema, The Electrical Double Layer Around a Spherical Colloid Particle, J.
Electrochem. Soc. 108 (12) (1961) 269C. doi:10.1149/1.2427992.
[89] D. Hänel, Molekulare Gasdynamik, Springer, 2004.
[90] S. Chapman, T. G. Cowling, The Mathematical Theory of Non-Uniform Gases, 3rd Edition, Cambridge Univ. Press,
1990.
[91] D. A. Wolf-Gladrow, Lattice-Gas Cellular Automata and Lattice Boltzmann Models: An Introduction, no. 1725 in
Lattice-gas Cellular Automata and Lattice Boltzmann Models: An Introduction, Springer, 2000.
[92] D. d’Humières, Generalized lattice-Boltzmann equations, in: Rarefied Gas Dynamics: Theory and Simulations, Vol. 159
of Prog. Astronaut. Aeronaut., AIAA, 1992, pp. 450–458. doi:10.2514/5.9781600866319.0450.0458.
[93] X. He, L.-S. Luo, Lattice Boltzmann Model for the Incompressible Navier-Stokes Equation, J. Stat. Phys. 88 (3) (1997)
927–944. doi:10.1023/B:JOSS.0000015179.12689.e4.
[94] Y. H. Qian, D. d’Humières, P. Lallemand, Lattice BGK Models for Navier-Stokes Equation, Europhys. Lett. 17 (6) (1992)
479. doi:10.1209/0295-5075/17/6/001.
[95] J. D. Sterling, S. Chen, Stability Analysis of Lattice Boltzmann Methods, J. Comput. Phys. 123 (1) (1996) 196–206.
doi:10.1006/jcph.1996.0016.
[96] L.-S. Luo, W. Liao, X. Chen, Y. Peng, W. Zhang, Numerics of the lattice Boltzmann method: Effects of collision models
on the lattice Boltzmann simulations, Phys. Rev. E 83 (5) (2011) 056710. doi:10.1103/PhysRevE.83.056710.
[97] P. L. Bhatnagar, E. P. Gross, M. Krook, A Model for Collision Processes in Gases. I. Small Amplitude Processes in
Charged and Neutral One-Component Systems, Phys. Rev. 94 (3) (1954) 511–525. doi:10.1103/PhysRev.94.511.

32

D. Bartuschat and U. Rüde / Journal of Computational Physics 00 (2017) 1–33

33

[98] M. Geier, M. Schönherr, A. Pasquali, M. Krafczyk, The cumulant lattice Boltzmann equation in three dimensions: Theory
and validation, Comput. Math. Appl. 70 (4) (2015) 507–547. doi:10.1016/j.camwa.2015.05.001.
[99] F. J. Higuera, S. Succi, R. Benzi, Lattice Gas Dynamics with Enhanced Collisions, Europhys. Lett. 9 (4) (1989) 345.
doi:10.1209/0295-5075/9/4/008.
[100] I. Ginzburg, F. Verhaeghe, D. d’Humières, Study of Simple Hydrodynamic Solutions with the Two-Relaxation-Times
Lattice Boltzmann Scheme, Commun. Comput. Phys. 3 (3) (2008) 519–581.
[101] L.-S. Luo, Analytic Solutions of Linearized Lattice Boltzmann Equation for Simple Flows, J. Stat. Phys. 88 (3/4) (1997)
913–926. doi:10.1023/B:JOSS.0000015178.19008.78.
[102] A. J. C. Ladd, R. Verberg, Lattice-Boltzmann Simulations of Particle-Fluid Suspensions, J. Stat. Phys. 104 (5) (2001)
1191–1251. doi:10.1023/A:1010414013942.
[103] Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E
65 (4) (2002) 046308. doi:10.1103/PhysRevE.65.046308.
[104] C. Rettinger, U. Rüde, A comparative study of fluid-particle coupling methods for fully resolved lattice Boltzmann
simulations, Comput. Fluids 154 (2017) 74–89. doi:10.1016/j.compfluid.2017.05.033.
[105] E. Fattahi, C. Waluga, B. Wohlmuth, U. Rüde, M. Manhart, R. Helmig, Lattice Boltzmann methods in porous media
simulations: From laminar to turbulent flow, Comput. Fluids 140 (2016) 247 – 259. doi:10.1016/j.compfluid.2016.
10.007.
[106] S. Bogner, S. Mohanty, U. Rüde, Drag correlation for dilute and moderately dense fluid-particle systems using the lattice
Boltzmann method, Int. J. Multiph. Flow 68 (2015) 71–79. doi:10.1016/j.ijmultiphaseflow.2014.10.001.
[107] P. Wesseling, Introduction To Multigrid Methods, R. T. Edwards, 2004.
[108] P. Knabner, L. Angermann, Numerical Methods for Elliptic and Parabolic Partial Differential Equations, Texts in Applied
Mathematics, Springer, 2003.
[109] B. I. des Poids et Mesures, The International System of Units (SI), 8th Edition, BIPM, 2006 [cited Oct. 2015].
URL http://www.bipm.org/en/publications/si-brochure/
[110] I. 80000-1:2009, Quantities and units – Part 1: General, ISO/IEC, 2009 [cited Oct. 2015].
URL https://www.iso.org/obp/ui/#iso:std:30669:en
[111] S. Donath, Wetting Models for a Parallel High-Performance Free Surface Lattice Boltzmann Method, Ph.D. thesis,
University of Erlangen-Nürnberg (2011).
URL http://www.dr.hut-verlag.de/978-3-8439-0066-9.html
[112] D. Bartuschat, D. Ritter, U. Rüde, Parallel multigrid for electrokinetic simulation in particle-fluid flows, in: High
Performance Computing and Simulation (HPCS). Madrid 2012, IEEE, 2012, pp. 374–380. doi:10.1109/HPCSim.2012.
6266940.
[113] R. Ramadugu, S. P. Thampi, R. Adhikari, S. Succi, S. Ansumali, Lattice differential operators for computational physics,
Europhys. Lett. 101 (5) (2013) 50006. doi:10.1209/0295-5075/101/50006.
[114] K. Iglberger, Software design of a massively parallel rigid body framework, Ph.D. thesis, University of Erlangen-Nürnberg
(2010).
URL http://www.dr.hut-verlag.de/978-3-86853-736-9.html
[115] L. Hufnagel, Transient simulation of electric double layers in electrokinetic flows with a coupled lattice Boltzmann
Link-Flux method, Bachelor thesis, Lehrstuhl für Informatik 10 (Systemsimulation), Universität Erlangen-Nürnberg
(November 2014).

33

