Mesoscale simulation with LB2D

Jonathan Chin <j.chin@qmul.ac.uk>

LB2D is a code for running lattice Boltzmann simulations of two-dimensional fluids, using the D2Q9 lattice. It supports porous media, two-component fluids with Shan-Chen interactions, and a variety of boundary conditions.

This document gives a brief description of the design of the LB2D code, followed by a tutorial on using it.

Overview of the code

Data structures

The majority of the data for the simulation are stored in a single simulation structure, defined in the code as the sim datatype. This is composed of several smaller components:

• A set of parameters which remain constant through a single simulation run, such as the lattice size and initial conditions.

• A set of parameters which may be varied ("steered") during the course of a simulation, such as the relaxation times of the fluids, forcing rates, or coupling constants.

• The lattice on which the simulation takes place.

The lattice structure is simply an array of site structures. Each site structure contains the entire state of a single site on the lattice: the red distribution function $f^r_i$ , the blue distribution function $f^b_i$ , and a flag to indicate whether or not there is an obstacle at the site. There is also a temporary copy of the lattice structure, which is used to hold data during the collision and advection timesteps.

Figure 1: Data structures in LB2D

Program flow

The core of the code is contained in the file lbe2d.c , which contains routines which operate on the sim structure, performing an advection or collision step, and imposing boundary conditions. Constants such as the lattice vectors are defined in lbe2d-const.c.

The file lbe2d-io.c contains routines to read input data (such as a description of a porous medium) from files, and to write simulation data to files.

The routines for reading and parsing an input file are contained in parse.c ; buffer.c , parse.c , and symbol.c are auxiliary libraries.

The overall program flow is handled in lbe2d-taskfarm.c. It reads and parses the input file, and allocates and initialises the simulation structure. The main loop of the program proceeds as follows:

• The forces (such as gravitational body force or Shan-Chen immiscibility force) are calculated and stored in a buffer in the lattice structure.

• A collision step is performed: at each site on the lattice, particles at that site are redistributed in a manner which conserves mass and momentum.

• A halo exchange is performed, imposing periodic boundary conditions.

• An advection step is performed, moving particles from one site to another. Particles are bounced off obstacle sites.

• Boundary conditions are imposed, such as ensuring constant fluid density or composition in certain regions.

• Another halo exchange is performed.

• Any required I/O takes place.

This repeats for a specified number of timesteps, or until an error condition occurs (such as if a negative particler density is detected). Once the simulation has finished, the memory for the lattice structure is released, and the program terminates.

Figure 2: Program flow in LB2D

Output file format

For very large production runs, a simulation will often be run on a large computer such as a Cray or SGI Origin, and its output will then be processed and examined on a smaller machine, such as a desktop workstation. This means that the simulation should produce output which can be read on any platform, not just the one on which it was run.

Unfortunately, most of the interesting scientific data from a simulation will be in floating-point format, and many computers have different and incompatible ways of storing floating-point data. Storing numbers as ASCII strings is an option, but wastes huge amounts of disk space and processor time.

To avoid this problem, LB2D writes its data in a format known as XDR ("eXternal Data Representation"), which was developed at Sun Microsystems as a standard way of exchanging data between different machines on a network. XDR forms the basis of the Network Filing System used to transfer files between UNIX machines. It has several advantages:

• XDR is an open standard, defined in RFC 1014. An RFC is an Internet Engineering Task Force document: other RFCs define the way data such as e-mail or web traffic are passed over the Internet.

• XDR is ubiquitous: libraries to convert data to and from the format can be found on almost any UNIX machine. Other portable formats such as NetCDF usually require the installation of separate libraries.

• XDR is very simple to use, and no royalties are payable for its use.

• XDR can store 64-bit binary floating-point numbers.

Scalar data

Scalar datasets, such as the density of a fluid component, or the order parameter during a phase transition, are written as arrays of double-precision XDR floating point numbers, with the X coordinate varying fastest.

A piece of the LB2D source code responsible for writing the order parameter (difference in red and blue densities) runs as follows:

        for (y=1;y<=NY;y++) {
                for (x=1;x<=NX;x++) {
        
                        phi=0.0;
        
                        for (s=0;s<nvecs;s++) {
                                phi+= SITE(x,y)->nr[s] - SITE(x,y)->nb[s];
                        }
        
                        xdr_double(&xdrs,&phi);
        
                }
        }

The lattice sites have X coordinates from 1 to NX -1, and Y coordinates from 1 to NY -1; coordinates 0 and NX or NY refer to the halo regions. SITE(x,y) is a macro which returns a pointer to the structure containing the state of the lattice site at coordinates $(x,y)$, and nvecs is the number of lattice vectors, set to 9 for the D2Q9 lattice. The innermost loop sums up the number of red particles at a site, and subtracts the number of blue particles; the final value is the order parameter which is written to a file by the xdr_double routine. The outer two loops move across the lattice.

Vector data

The velocity field of a fluid is a two-dimensional vector. The velocity output files contain pairs of XDR-double precision floating point numbers, written in the same order as the scalar data. The first element in each pair is the X velocity; the second element is the Y velocity.

Running LB2D

A single simulation can be run by constructing an input file called input-0000 , and running the lbe2d-taskfarm program. On parallel machines, the program can be made to run several simultaneous independent simulations (a "task farm") from consecutively-numbered input files.

The simulations will typically write some data such as a velocity or density field to files; these files are then examined once the simulation has completed.

Input file example

An example of a valid input file is:

        { folder="poiseuille"; makedir=1;
        
            initcond="allblue";
            wallsides=1;
            writeyvel=1;
            tau_b = 0.547 ; rho = 1.0;
            tmax=100000; dt=10000;
            gravity=0.00001;
            nx=32 ; ny=64;
        
        }

This file produces a simulation of single-component Poiseuille flow through a channel, which runs for 100000 timesteps, and writes the Y-velocity of the fluid every 10000 timesteps.

A single simulation is described by a set of key-value pairs enclosed in braces. Each pair must end with a semicolon, and whitespace is ignored. String variables must be enclosed in double quotes. Boolean variables should be set to 1 to indicate true, and 0 to indicate false.

folder="poiseuille";

This instructs LB2D to write all its output files into a directory called poiseuille.

initcond="allblue";

Sets the initial state of the system to contain only the blue component, at rest, at a density given by the rho parameter (see below).

makedir=1;

Instructs LB2D to create the poiseuille directory, specified above, if it does not already exist.

wallsides=1;

Creates obstacle walls on the left and right sides (low X and high X) of the lattice, producing a two-dimensional channel.

writeyvel=1;

This Boolean parameter turns on output of the Y-velocity during the I/O stage.

tau_b = 0.547;

Sets the relaxation time $\tau_b$ of the blue component to $0.547$, giving a viscosity of $(\tau-\frac 1 2 )/3 \simeq 0.015667$ in lattice units.

rho = 1.0;

Sets the initial density of blue particles to $1.0$ particles per site.

tmax=100000;

Terminates the simulation after $100000$ timesteps.

dt=10000;

Sets the output to be written every $10000$ timesteps.

gravity=0.00001;

Sets the body force acting on all components (i.e.the gravitational acceleration) to take the value $0.00001$ in lattice units.

nx=32 ; ny=64;

Sets the width and height of the simulation lattice.

List of possible input file parameters

I/O and global parameters

nx (integer)

Sets the X size of the simulation lattice.

ny (integer)

Sets the Y size of the simulation lattice.

tmax (integer)

Sets the number of timesteps for which the simulation will run.

dt (integer)

Sets the number of timesteps between which the code will write output files.

folder (STRING)

The name of a directory into which the code will move before running a simulation. If running many simultaneous simulations, this is useful for keeping the outputs from each separate run in separate directories. If the makedir option is set, then this directory will be created if it does not already exist.

makedir (Boolean)

If set to 1, then the directory specified by the folder option will be created if it does not already exist.

writecolour (Boolean)

If set, the colour field will be written to a new file every dt timesteps.

writered (Boolean)

If set, the red particle density will be written to a new file every dt timesteps.

writeblu (Boolean)

If set, the blue particle density will be written to a new file every dt timesteps.

writeyvel (Boolean)

If set, the Y velocity will be written to a new file every dt timesteps.

writept (Boolean)

If set, the pressure tensor will be written to a new file every dt timesteps.

Simulation parameters

tau_r,tau_b (float)

Set the relaxation time of the red component and blue component, respectively.

g_cc (float)

Sets the value of the Shan-Chen coupling constant between the red and blue components.

gravity (float)

Sets the value of a body force acting in the Y direction.

m_r,m_b (float)

Sets the value of the mass of the red and blue particles.

Initial condition parameters

initcond (STRING)

Sets the initial conditions for the simulation. Can take one of the following values:

bubble - creates a bubble of red component in a sea of blue. The radius of the bubble is specified in the radius option; the particle densities in the system are specified by the rho option.

cone - like bubble, but the inside of the bubble is not pure red; instead, the proportion of red is linearly increased and the proportion of blue linearly decreased closer to the centre of the bubble.

capwave - creates a capillary wave between regions of homogeneous red and blue component. The amplitude and wavelength of the wave are given by the amp and lambda options; the density of particles is given by the rho option.

flatface - like capwave, but produces a straight-line interface.

channel - assumes there are bounceback walls at x=1 and x=NX, and produces a channel of red component inside regions of blue component. If the system has an even number of points in the X direction, then the interface will be on a lattice site; in this case, the interfacial sites are given a 50/50 mixture.

allblue - fills the simulation arena with blue component.

allred - fills the simulation arena with red component.

finger - sets up a viscous fingering simulation: the system is filled with blue component with an equilibrium Poiseuille velocity profile as calculated from the gravity option.

spinodal - sets up a spinodal decomposition simulation. The system is filled with a 50/50 mixture of red and blue, with the total density given by the rho option. This is then perturbed by an amount given by the perturb option to avoid metastable initial states.

rockfile (string)

If set, specifies the name of a file containing rock geometry data, which is loaded after the non-rock-related initial conditions have been set. If the rockfileformat option is not set, then it is assumed that the rock file is a greyscale PNG image.

rockfileformat (string)

Sets the format of the rock geometry file named in the rockfile option. Possible values are "xdr", for XDR format, "pgm" for a Portable GreyMap image, or "png" for a Portable Network Graphics (PNG) image. It is recommended to use PNG, since this format is widely supported by image programs, and uses lossless compression to reduce the size of the files.

wallsides (Boolean)

If set, produces rock walls on the low X and high X extremities of the simulation region, to facilitate simulation of channel flow. This happens after any rock data from the rockxdr option is loaded.

rho (float)

Sets the initial particle density.

perturb (float)

Once the other initial conditions have been set, this will perturb the density of rest particles at each site on the lattice by the given fraction. For example, if the lattice is initialised to a constant density of 1.0 at each site, then setting perturb=0.1 will give a lattice with particles densities ranging from 0.9 to 1.1.

radius (float)

For the bubble initial condition, the bubble is given a radius of radius*(smallest system dimension). For example, to produce a bubble of radius 64 on a 256x256 lattice, set radius=0.25.

lambda (float)

For capillary wave initial conditions, sets the wavelength in lattice sites of the initial waveform.

rockr,rockb (float)

Set the red and blue particle density at rock sites. If rockr=rockb=0.0, then rock sites are neutrally wetting.

amp (float)

For capillary wave initial conditions, sets the amplitude in lattice sites of the initial waveform.

seed (integer)

Sets the random number generator seed. For ensemble averaging, this should be set to a different value for each run un the ensemble.

Example simulations

Poiseuille flow

The Navier-Stokes equation for fluid flow is:

$$\mathbf{\dot{u}} + (\mathbf{u} \cdot \mathbf{\nabla} )\mathbf{u} + \frac{1}{\rho} \mathbf{\nabla} p = \nu \nabla^2\mathbf{u}$$

$\mathbf{u}$ is the fluid velocity vector; $\rho$ the fluid density, $p$ the scalar pressure, and $\nu$ the kinematic viscosity. Consider a 2D channel, bounded on either side by no-slip walls. The pressure gradient $\nabla p$ is equal to the applied body force $-G$ . Assuming that $d / d y = 0$ and $\mathbf{u} = (u_x,u_y)=(0,u)$ simplifies the equation to give

$$\frac{d^2 u}{d x^2} = - \frac{G}{\rho \nu}$$

This solves to give a quadratic velocity profile:

$$u = \frac{G}{2 \rho \nu} \left( \frac{L^2}{4} - x^2 \right)$$

Setting up a simulation

An input file to produce Poiseuille flow is:

        {
        
            folder="poiseuille"; makedir=1;
            initcond="allblue";
            rho=1.0;
            wallsides=1;
            writeyvel=1;
            tau_b=0.547; gravity=0.00001;
            tmax=10000; dt=1000;
            nx=32; ny=64;
        
        }

This should be entered (using Notepad or any other text editor) into a plain ASCII file called input-0000, in the same directory as the lbe2d-forker executable.

Running a simulation

Once the input file has been set up, open a Cygwin shell, and tell LB2D to launch a single simulation:

         % ./lbe2d-forker 1

All being well, the code should print something like the following:

lbe2d-visco taskfarm version
$Id: lbetute.tex,v 1.9 2002/07/11 14:06:56 jon Exp $
This is rank 0.
Initializing to all blue, rho 1.000000.
Creating side walls.

{
folder = "poiseuille";
makedir = 1;
initcond = "allblue";
rho = 1.000000;
wallsides = 1;
writeyvel = 1;
tau_b = 0.547000;
gravity = 0.000010;
tmax = 10000;
dt = 1000;
nx = 32;
ny = 64;
}

Running...
Red: 0.0000000000 Blue: 1920.0000000001 Sum: 1920.0000000001 px: 0.0000000000
py: 0.0187733333
Timestep 0
Timestep 100
Timestep 200
Timestep 300
Timestep 400

Every time output is written to a file, the total density and momentum of particles is also written to the screen. In the absence of obstacles, the total momentum (px and py) should be a constant; in this case, the body force adds momentum to the system, which is then taken away by collisions with the side walls. Hence, the X momentum is constant and the Y momentum increases, settling down to a constant value once the velocity profile has stabilised.

Examining simulation output

Once the simulation has finished, the poiseuille directory should be populated with files. Files with the extension .xdr contain numerical data; files with the extension .fld describe the size and format of the XDR files so that they may be imported into other applications.

To examine the velocity profile at $x = 32$ at timestep 10000, change into the poiseuille directory and type:

        % ../xdrslice 32 32 yvel*10000.xdr 

This should print the velocity profile to the screen; this can be redirected to a file and fed into a graph plotting program such as Gnuplot. The development of the quadratic profile, and comparison with the theoretical form of the profile are shown below.

Figure 3: Evolution of the velocity profile towards the quadratic Poiseuille profile

Figure 4: Poiseuille velocity profile.

Spinodal decomposition in 2D

If two fluids are mixed, and then the temperature is suddenly reduced, so that they become immiscible, they separate in a process known as spinodal decomposition. This can be simulated in an LB model: the red and blue particles are made immiscible through the Shan-Chen interaction, and the initial condition is set to be a 50-50 mixture of the two phases. Note that there must be some slight fluctuations in the initial conditions, since a completely homogeneous 50-50 mixture would be metastable, and would not separate until a slight fluctuation is introduced.

Such an initial condition is built in to the LB2D code, and can be selected by setting the initcond variable to the value "spinodal", as in the following input file:

        { folder="spinodal";
        initcond="spinodal";
        writecolour="true";
        makedir = "yes";
        g_cc=1.08 ; tau_r = 1.0 ;tau_b = 1.0; rho = 1.0;
        tmax=5000 ; dt = 50 ; gravity=0.0;
        nx=256 ; ny=256;
        }

The structure factor can be calculated using the 2dsf program: to calculate the averaged structure factor at timestep 1000, type:

        % ../2dsf 256 256 col*t01000.xdr

Examination of the structure factor at successive timesteps should show a curve whose overall shape remains the same, but grows exponentially in amplitude.

Figure 5: Poiseuille velocity profile.

Modifying the collision operator

Implementation of the collision operator

At each site, the collision process has five main steps:

• Calculation of the density $\rho^{\sigma} = \sum_{\sigma} f^{\sigma}_i$ and momentum $\mathbf{u}^{\sigma} = \sum_{\sigma} f_i^{\sigma}\mathbf{c}_i$ for each component $\sigma$. For a two-component model, $\sigma$ refers to either red or blue.

• Calculation of the weighted average velocity $\mathbf{u} = \left( \sum _{\sigma} \rho^{\sigma}\mathbf{u}^{\sigma} / \tau^{\sigma} \right) / \left( \sum_{\sigma} \rho^{\sigma} / \tau^{\sigma} \right)$

• Calculation of the perturbed velocity $\mathbf{v}^{\sigma} =\mathbf{u} + \tau^{\sigma}\mathbf{F}^{\sigma} / \rho^{\sigma}$ from the force $\mathbf{F}^{\sigma}$ acting on each component.

• Calculation of the equilibrium distribution function:
  

  $$f_i^{\mathrm{eq}\; \sigma} = w_i n^{\sigma} \left[1 + 3 \l... ...ight) - \frac{3}{2} \vert\mathbf{v}^{\sigma} \vert^2 \right]$$
  
  

• Collision via the BGK operator:

  
  

  $$f^{' \sigma}_i = f^{\sigma}_i - \frac{1}{\tau} \left( f^{\sigma}_i - f^{\mathrm{eq}\; \sigma}_i \right)$$
  
  

This is all performed in the routine sim_collide in the file lbe2d.c . This routine is called with a pointer to the simulation data structure in the variable mySim.

Simple floating-point simulation parameters, such as relaxation times, can be found by examining individual components of this structure: for example, the relaxation time of the blue component is given my the expression mySim->tau_b ; the Shan-Chen coupling constant, used to determine the strength of the force between two components, is given by mySim->g_cc .

To facilitate working with the code, several preprocessor macros are defined to extract further information from the simulation data structure:

• SITE(x,y) returns a pointer to the lattice site at coordinates $(x,y)$. SITE(x,y)->nr[i] gives the red-component distribution function for the ith lattice vector at site $(x,y)$, i.e. $f_i^r ( x , y )$. Similarly, SITE(x,y)->nb[i] gives the blue distribution function, and SITE(x,y)->rock returns the value 1 if there is an obstacle at lattice site $(x,y)$.

• FORCE(x,y) describes the external forces acting on the fluid at $(x,y)$. The $x$ component of the force on the red component at $(x,y)$ is given by FORCE(x,y)->red.x , the $y$ component by FORCE(x,y)->red.y , and similarly for FORCE(x,y)->blu .

The collision code begins with a loop over all lattice sites:

        for (y=1;y{<}=NY;y++) {
        
                for (x=1;x{<}=NX\textrm{};x++) {

[...]

                }
        }

NX and NY are preprocessor macros which determine the dimensions of the lattice from the data in the mySim structure. For brevity and simplicity of code, the value of SITE(x,y) is stored in a variable called me inside the loop, so that, for example, me->nr[i] gives the red density on the $i$th lattice vector of the current site. The new values of the distribution function cannot be stored straight back into the lattice, since this would lead to incorrect calculation of the distribution function at other stites. Instead, the new values are stored in a temporary array, the lattice sites of which can be accessed through the TMPSITE(x,y) macro. Inside the loop, the current value of TMPSITE(x,y) is stored in a variable called him.

The first stage of collision calculates the velocities and densities by summing over all the lattice vectors:

        nr = nb = 0.0;
        
        upx_r = upy_r = upx_b = upy_b = 0.0;
        
        for (s=0;s<nvecs;s++) {
        
                nr += me->nr[s];
                nb += me->nb[s];
                
                upx_r += me->nr[s]*e[s][0];
                upy_r += me->nr[s]*e[s][1];
                upx_b += me->nb[s]*e[s][0];
                upy_b += me->nb[s]*e[s][1];
        
        }
        
        rho_r = nr * m_r;
        rho_b = nb * m_b;

Now, the red component momentum is in upx_r and upy_r and the red density is in rho_r, similarly for the blue component. Then, the equilibrium velocity is calculated:

        ux_r = ( m_r/tau_r * upx_r + m_b/tau_b * upx_b) / 
                ( rho_r/tau_r + rho_b/tau_b );

        uy_r = ( m_r/tau_r * upy_r + m_b/tau_b * upy_b) / 
                ( rho_r/tau_r + rho_b/tau_b );

This has the same value for each component, so it is only calculated once, and copied into the ux_b and uy_b variables. The next step is to perturb this velocity according to the Shan-Chen forcing term, which has already been calculated:

        Fx_r = FORCE(x,y)->red.x;
        Fy_r = FORCE(x,y)->red.y;

        if (0.0 < rho_r) {
                ux_r += tau_r/rho_r * Fx_r;
                uy_r += tau_r/rho_r * Fy_r;
        }

Note that this only occurs if the component has nonzero density: in a region containing only blue component, rho_r would be zero, resulting in a divide-by-zero error in an otherwise physically sane situation.

Before calculating the equilibrium distribution, some frequently-used values are precalculated: while modern compilers should efficiently deal with common subexpressions, this technique will often make the code easier to read and work with. For the red component, the dot product $\mathbf{v}^{\sigma} \cdot \mathbf{c}\_i$ is stored in edotu_r , and the squared modulus of $\mathbf{v}^{\sigma}$ in u2_r . The equilibrium distribution is then straightforward:

        feq_r = w[a] * nr * ( 
                1.0 + 
                3.0*edotu_r + 4.5*edotu_r*edotu_r -
                1.5*u2_r );

The values w[i] are the weights required to produce isotropic hydrodynamics on a D2Q9 lattice, and are stored in a constant array defined at compile-time. Once the equilibrium distribution is known, all that remains is for the BGK collision to occur:

        him->nr[a] = (1.0-1.0/tau_r)*(me->nr[a]) + 1.0/tau_r*feq_r;

Modification of the collision operator

During the process of phase separation of a mixture of fluids, hydrodynamic processes play an important role: instead of the separation being governed by the random diffusion of one fluid through another, directed flows of a component from smaller domains towards larger ones may occur.

Hydrodynamics may be crudely be disabled by setting the equilibrium velocity $\mathbf{u}$ to zero during the collision step. This means that the only force acting on the components is the Shan-Chen force, and that momentum cannot be transferred between lattice sites during the advection step.

Load

lbe2d.c

into a text editor, and find the lines containing the code to calculate the value of $\mathbf{u}$ :

        ux_r = ( m_r/tau_r * upx_r + m_b/tau_b * upx_b) / 
                ( rho_r/tau_r + rho_b/tau_b );
        uy_r = ( m_r/tau_r * upy_r + m_b/tau_b * upy_b) / 
                ( rho_r/tau_r + rho_b/tau_b );

Delete these, and replace with the line:

        ux_r=uy_r=0.0;

Then, save the file, and run make to rebuild the simulation program. If this is successful, then re-run the spinodal decomposition simulation, making sure to change the folder option so as not to write over previous results. The components should separate much slower, with long and thin domains gradually thickening and coalescing, in comparison to the hydrodynamic regimes which quickly form large droplets.

About this document ...

Mesoscale simulation with LB2D

This document was generated using the LaTeX2HTML translator Version 2002 (1.1)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -split 0 lbetute

The translation was initiated by Jonathan Chin on 2003-07-14