With this tutorial you can get started using the Lattice-Boltzmann method
for scientific applications. We give a brief introduction about the theory
and how to use it in ESPResSo. We have selected three interesting problems for which LB can be applied and which are well understood. You can start with any of them.
The tutorial is relatively long and working through it carefully is work for at least a full day. You can however get a glimpse of different aspects by starting to work on the tasks.
Note: LB can not be used as a black box. It is unavoidable to spend time
learning the theory and gaining practical experience.
In this tutorial, you will learn basics about the Lattice Boltzmann Method (LBM) with
special focus on the application on soft matter simulations or more precisely on how to
apply it in combination with molecular dynamics to take into account hydrodynamic
solvent effects without the need to introduce thousands of solvent particles.
The LBM – its theory as well as its applications – is still a very active field of research.
After almost 20 years of development there are many cases in which the LBM has proven
to be fruitful, in other cases the LBM is considered promising, and in some cases it has
not been of any help. We encourage you to contribute to the scientific discussion of
the LBM because there is still a lot that is unknown or only vaguely known about this
fascinating method.
This tutorial should enable you to start a scientific project applying the LB method with ESPResSo. In the first part we summarize a few basic ideas behind LB and describe the interface. In the second part we suggest three different classic examples where hydrodynamics are important. These are
Hydrodynamic resistance of settling particles. We measure the drag force of single particles and arrays of particles when sedimenting in solution.
Polymer diffusion. We show that the diffusion of polymers is accelerated by hydrodynamic interactions.
Poiseuille flow. We reproduce the flow profile between two walls.
With Version 3.1 ESPResSo has learned GPU support for LB. We recommend however
version 3.3, to have all features available. We absolutely recommend using the GPU code,
as it is much (100x) faster than the CPU code.
For the tutorial you will have to compile in the following features: PARTIAL_PERIODIC,
EXTERNAL_FORCES, ELECTROSTATICS, LB_GPU, LB_BOUNDARIES_GPU,
LENNARD_JONES.
All necessary files for this tutorial are located in the directory
Here we want to repeat a few very basic facts about the LBM. You will find much better
introductions in various books and articles, e.g. [1, 2]. It will however help clarifying
our choice of words and we will eventually say something about the implementation in
ESPResSo. It is very loosely written, with the goal that the reader understands basic
concepts and how they are implemented in ESPResSo.
The LBM essentially consists of solving a fully discretized version of the linearized
Boltzmann equation. The Boltzmann equation describes the time evolution of the one
particle distribution function $f (x, p, t)$, which is the probability to find a molecule in
a phase space volume $dxdp$ at time $t$.The function $f$ is normalized so that the integral
over the whole phase space is the total mass of the particles:
$$\int f (x, p) dxdp = N m,$$
where $N$ denotes the particle number and m the particle mass. The quantity $f(x, p) dxdp$
corresponds to the mass of particles in this particular cell of the phase space, the
population.
The LBM discretizes the Boltzmann equation not only in real space (the lattice!) and
time, but also the velocity space is discretized. A surprisingly small number of velocities,
in 3D usually 19, is sufficient to describe incompressible, viscous flow correctly. Mostly
we will refer to the three-dimensional model with a discrete set of 19 velocities, which is
conventionally called D3Q19. These velocities, $\vec{c_i}$ , are chosen so that they correspond to
the movement from one lattice node to another in one time step. A two step scheme is
used to transport information through the system: In the streaming step the particles
(in terms of populations) are transported to the cell where they corresponding velocity
points to. In the collision step, the distribution functions in each cell are relaxed towards
the local thermodynamic equilibrium. This will be described in more detail below.
The hydrodynamic fields, the density, the fluid momentum density, the pressure tensor can be calculated straightforwardly from the populations. They correspond to the
moments of the distribution function:
Here the Greek indices denotes the cartesian axis and the Latin indices indicate the number in the disrete velocity set. Note that the pressure tensor is symmetric. It is easy to see that these equations are linear transformations of the $f_i$ and that they carry the most important information. They are 10 independent variables, but this is not enough to store the full information of 19 populations. Therefore 9 additional quantities are introduced. Together they form a different basis set of the 19-dimensional population space, the modes space and the modes are denoted by $m_i$. The 9 extra modes are referred to as kinetic modes or ghost modes. It is possible to explicitly write down the base transformation matrix, and its inverse and in the ESPResSo LBM implementation this basis transformation is made for every cell in every LBM step. It is possible to write a code that does not need this basis transformation, but it has been shown, that this only costs 20% of the computational time and allows for larger flexibility.
The second step is the collision part, where the actual physics happens. For the LBM it is assumed that the collision process linearly relaxes the populations to the local equilibrium, thus that it is a linear (=matrix) operator acting on the populations in each LB cell. It should conserve the particle number and the momentum. At this point it is clear why the mode space is helpful. A 19 dimensional matrix that conserves the first 4 modes (with the eigenvalue 1) is diagonal in the first four rows and columns. Some struggling with lattice symmetries shows that four independent variables are enough to characterize the linear relaxation process so that all symmetries of the lattice are obeyed. Two of them are closely related to the shear and bulk viscosity of the fluid, and two of them do not have a direct physical equivalent. They are just called relaxation rates of the kinetic modes.
The equilibrium distribution to which the populations relax is obtained from maximizing the information entropy $\sum f_i \log f_i$ under the constraint that the density and velocity take their particular instantaneous values.
In mode space the equilbrium distribution is calculated much from the local density and velocity. The kinetic modes 11-19 have the value 0 in equilibrium. The collision operator is diagonal in mode space and has the form \begin{align*} m^\star_i &= \gamma_i \left( m_i - m_i^\text{eq} \right) + m_i ^\text{eq}. \end{align*} Here $m^\star_i$ is the $i$th mode after the collision. In words we would say; each mode is relaxed towards it's equilibrium value with a relaxation rate $\gamma_i$. The conserved modes are not relaxed, or, the corresponding relaxation parameter is one.
By symmetry consideration one finds that only four independent relaxation rates are allowed. We summarize them here. \begin{align*} m^\star_i &= \gamma_i m_i \\ \gamma_1=\dots=\gamma_4&=1 \\ \gamma_5&=\gamma_\text{b} \\ \gamma_6=\dots=\gamma_{10}&=\gamma_\text{s} \\ \gamma_{11}=\dots=\gamma_{16}&=\gamma_\text{odd} \\ \gamma_{17}=\dots = \gamma_{19}&=\gamma_\text{even} \\ \end{align*}
To include hydrodynamic fluctuations of the fluid, random fluctuations are added to the nonconserved modes $4\dots 19$ on every LB node so that the LB fluid temperature is well defined and the corresponding fluctuation formula, according to the fluctuation dissipation theorem holds. An extensive discussion of this topic is found in [3].
Particles are coupled to the LB fluid with the force coupling; the fluid velocity at the position of a particle is calculated by a multilinear interpolation and a force is applied on the particle that is proportional to the velocity difference between particle and fluid: \begin{equation} \vec{F}_D = - \gamma \left(v-u\right) \end{equation} The opposite force is distributed on the surrounding LB nodes. Additionally a random force is added to maintain a constant temperature, according to the fluctuation dissipation theorem.
ESPResSo features two virtually independent implementations of LB. One implementation uses CPUs and one uses a GPU to perform the computational work. For this, we provide two actor classes LBFluid and LBFluid_GPU with the module espressomd.lb,
as well as the optional LBBoundary class found in espressomd.lbboundaries.
The LB lattice is a cubic lattice, with a lattice constant agrid that
is the same in all spacial directions. The chosen box length must be an integer multiple
of agrid. The LB lattice is shifted by 0.5 agrid in all directions: the node
with integer coordinates $\left(0,0,0\right)$ is located at
$\left(0.5a,0.5a,0.5a\right)$.
The LB scheme and the MD scheme are not synchronized: In one
LB time step, several MD steps may be performed. This allows to speed
up the simulations and is adjusted with the parameter tau.
The LB parameter tau must be an integer multiple of the MD timestep.
Even if MD features are not used, the System parameters cell_system.skin and time_step must be set. LB steps are performed
in regular intervals, such that the timestep $\tau$ for LB is recovered.
Important Notice: All commands of the LB interface use
MD units. This is convenient, as e.g. a particular
viscosity can be set and the LB time step can be changed without
altering the viscosity. On the other hand this is a source
of a plethora of mistakes: The LBM is only reliable in a certain
range of parameters (in LB units) and the unit conversion
may take some of them far out of this range. So note that you always
have to assure that you are not messing with that!
One brief example: a certain velocity may be 10 in MD units.
If the LB time step is 0.1 in MD units, and the lattice constant
is 1, then it corresponds to a velocity of $1\ \frac{a}{\tau}$ in LB units.
This is the maximum velocity of the discrete velocity set and therefore
causes numerical instabilities like negative populations.
The LBFluid class provides an interface to the LB-Method in the ESPResSo core. When initializing an object, one can pass the aforementioned parameters as keyword arguments. Parameters are given in MD units. The available keyword arguments are:
dens : The density of the fluid.
agrid : The lattice constant of the fluid. It is used to determine the number of LB nodes per direction from box_l. They have to be compatible.
visc : The kinematic viscosity
tau : The time step of LB. It has to be an integer multiple of System.time_step.
fric : The frction coefficient $\gamma$ for the coupling scheme.
ext_force : An external force applied to every node. This is given as a list, tuple or array with three components.
Using these arguments, one can initialize an LBFluid object. This object then needs to be added to the system's actor list. The code below provides a minimum example.
from espressomd import lb, assert_features
assert_features(['LB_GPU', 'LB_BOUNDARIES_GPU'])
# initialize the System and set the necessary MD parameters for LB.
system = System()
system.box_l = [31, 41, 59]
system.time_step = 0.01
system.cell_system.skin = 0.4
# Initialize and LBFluid with the minimum set of valid parameters.
lbf = lb.LBFluid_GPU(agrid = 1, dens = 10, visc = .1, tau = 0.01)
# Activate the LB by adding it to the System's actor list.
system.actors.add(lbf)The LBFluid class also provides a set of methods which can be used to sample data from the fluid nodes. For example
lbf[X ,Y ,Z].quantityreturns the quantity of the node with $(X, Y, Z)$ coordinates. Note that the indexing in every direction starts with 0. The possible properties are:
velocity : the fluid velocity (list of three floats)
pi : the stress tensor (list of six floats: $\Pi_{xx}$, $\Pi_{xy}$, $\Pi_{yy}$, $\Pi_{xz}$, $\Pi_{yz}$, and $\Pi_{zz}$)
pi_neq : the nonequilibrium part of the stress tensor, components as above.
population : the 19 populations of the D3Q19 lattice.
boundary : the boundary index.
density : the local density.
The LBBoundary class represents a boundary on the LBFluid
lattice. It is dependent on the classes of the module espressomd.shapes as it derives its geometry from them. For the initialization, the arguments shape and velocity are supported. The shape argument takes an object from the shapes module and the velocity argument expects a list, tuple or array containing 3 floats. Setting the velocity will results in a slip boundary condition.
Note that the boundaries are not constructed through the periodic boundary. If, for example, one would set a sphere with its center in one of the corner of the boxes, a sphere fragment will be generated. To avoid this, make sure the sphere, or any other boundary, fits inside the original box.
This part of the LB implementation is still experimental, so please tell us about your experience with it. In general even the simple case of no-slip
boundary is still an important research topic in the lb community and in combination with point particle coupling not much experience exists. This means: Do research on that topic, play around with parameters and figure out what happens.
Boundaries are initialized by passing a shape object to the LBBoundary class. One way to initialize a wall is:
from espressomd import lbboundaries, shapes
wall = lbboundaries.LBBoundary(shape=shapes.Wall(normal = [1, 0, 0], dist = 1), velocity = [0, 0, 0.01])
system.lbboundaries.add(wall)Note that all used variables are inherited from previous examples. This will create a wall a surface normal of $(1, 0, 0)$ at a distance of 1 from the origin of the coordinate system in direction of the normal vector. The wall exhibits a slip boundary condition with a velocity of $(0, 0, 0.01)$. For the a no-slip condition, leave out the velocity argument or set it so zero. Please refer to the user guide for a complete list of constraints.
Currently only the so called link bounce back method is implemented, where the effective hydrodynamic boundary is located midway between two nodes. This is the simplest and yet a rather effective approach for boundary implementation. Using the shape objects distance function, the nodes determine once during initialisation whether they are boundary or fluid nodes.
As a first test, we measure the drag force on different objects in a simulation box. Under low Reynolds number conditions, an object with velocity $\vec{v}$ experiences a drag force $\vec{F}_\text{D}$ proportional to the velocity: \begin{align*} \vec{F}_\text{D}=-\gamma\vec{v}, \end{align*} where $\gamma$ is denoted the friction coefficient. In general $\gamma$ is a tensor thus the drag force is generally not parallel to the velocity. For spherical particles the drag force is given by Stokes' law: \begin{align*} \vec{F}_\text{D}=-6\pi\eta a\vec{v}, \end{align*} where $a$ is the radius of the sphere.
In this task you will measure the drag force on falling objects with LB and ESPResSo. In the sample script lb_stokes_force.py a spherical object at rest is centered in a square channel.
In [ ]:
#This code is from 'lb_stokes_force.py'.
from espressomd import System, lb, lbboundaries, shapes
import numpy as np
import sys
# System setup
system = System(box_l = [1.0, 1.0, 1.0]) # The box_l is specified below.
# [1.0,1.0,1.0] is dummy values
agrid = 1
radius = 5.5
box_width = 64
real_width = box_width+2*agrid
box_length = 64
system.box_l = [real_width, real_width, box_length]
system.set_random_state_PRNG()
np.random.seed(seed = system.seed)
system.time_step = 0.2
system.cell_system.skin = 0.4
# The temperature is zero.
system.thermostat.set_lb(kT=0)
Bounce back boundary conditions are assumed on the sphere. At the channel boundary the velocity is fixed by using appropriate boundary conditions. Within a few hundred or thousand integration steps a steady state develops and the force on the sphere converges.
In [ ]:
# LB Parameters
v = [0,0,0.01] # The boundary slip
kinematic_visc = 1.0
# Invoke LB fluid
lbf = lb.LBFluidGPU(visc=kinematic_visc, dens=1
, agrid=agrid, tau=system.time_step, fric=1)
system.actors.add(lbf)
# Setup walls
walls = [None] * 4
walls[0] = lbboundaries.LBBoundary(shape=shapes.Wall(normal=[-1,0,0],
dist = -(1+box_width)), velocity=v)
walls[1] = lbboundaries.LBBoundary(shape=shapes.Wall(normal=[1,0,0],
dist = 1),velocity=v)
walls[2] = lbboundaries.LBBoundary(shape=shapes.Wall(normal=[0,-1,0],
dist = -(1+box_width)), velocity=v)
walls[3] = lbboundaries.LBBoundary(shape=shapes.Wall(normal=[0,1,0],
dist = 1),velocity=v)
for wall in walls:
system.lbboundaries.add(wall)
# setup sphere without slip in the middle
sphere = lbboundaries.LBBoundary(shape=shapes.Sphere(radius=radius,
center = [real_width/2] * 2 + [box_length/2],
direction = 1))
system.lbboundaries.add(sphere)
lbf.print_vtk_boundary("./boundary.vtk")
def size(vector):
tmp = 0
for k in vector:
tmp+=k*k
return np.sqrt(tmp)
for i in range(10):
system.integrator.run(400)
lbf.print_vtk_velocity("fluid%04i.vtk" %i)
sys.stdout.write("\rloop: %02i"%i)
sys.stdout.flush()
print("\nIntegration finished.")
# get force that is exerted on the sphere
force = sphere.get_force()
print("Measured force: f=%f" %size(force))
stokes_force = 6*np.pi*kinematic_visc*radius*size(v)
print("Stokes' Law says: f=%f" %stokes_force)
Measure the drag force for three different input radii of the sphere. How good is the agreement with Stokes’ law? Calculate an effective radius from Stokes’ law and the drag force measured in the simulation. Is there a clear relation to the input radius? Remember how the bounce back boundary condition work and how good spheres can be represented by them.
The script produces vtk files of the flow field. Visualize the flow field with paraview. Open paraview by typing it on the command line. Make sure you are in the folder where the files are located. So the agenda is:
Measure the drag force for a fixed radius but varying system size. Does the drag force increase of decrease with the system size? Can you find a qualitative explanation?
[1] S Succi. The lattice Boltzmann equation for fluid dynamics and beyond. Clarendon Press, Oxford, 2001.
[2] B. Dünweg and A. J. C. Ladd. Advanced Computer Simulation Approaches for Soft Matter Sciences III, chapter II, pages 89–166. Springer, 2009.
[3] B. Dünweg, U. Schiller, and A.J.C. Ladd. Statistical mechanics of the fluctuating lattice-boltzmann equation. Phys. Rev. E, 76:36704, 2007.
[4] P. G. de Gennes. Scaling Concepts in Polymer Physics. Cornell University Press, Ithaca, NY, 1979.
[5] M. Doi. Introduction do Polymer Physics. Clarendon Press, Oxford, 1996.
[6] Michael Rubinstein and Ralph H. Colby. Polymer Physics. Oxford University Press, Oxford, UK, 2003.
[7] Daan Frenkel and Berend Smit. Understanding Molecular Simulation. Academic Press, San Diego, second edition, 2002.
In [ ]:
In [ ]: