LBM Examples
2D flow around a cylinder
This example introduces the LBM module in 2D for a flow around a cylinder in 2d. The example code is as follows:
struct UserData
{
Array<Cell *> Left;
Array<Cell *> Right;
Array<double> Vel;
double vmax;
double rho;
};
void Setup (LBM::Domain & dom, void * UD)
{
UserData & dat = (*static_cast<UserData *>(UD));
// Cells with prescribed velocity
#ifdef USE_OMP
#pragma omp parallel for schedule(static) num_threads(dom.Nproc)
#endif
for (size_t i=0; i<dat.Left.Size(); ++i)
{
Cell * c = dat.Left[i];
if (c->IsSolid) continue;
double rho = (c->F[0]+c->F[2]+c->F[4] + 2.0*(c->F[3]+c->F[6]+c->F[7]))/(1.0-dat.Vel[i]);
c->F[1] = c->F[3] + (2.0/3.0)*rho*dat.Vel[i];
c->F[5] = c->F[7] + (1.0/6.0)*rho*dat.Vel[i] - 0.5*(c->F[2]-c->F[4]);
c->F[8] = c->F[6] + (1.0/6.0)*rho*dat.Vel[i] + 0.5*(c->F[2]-c->F[4]);
c->Rho = c->VelDen(c->Vel);
}
// Cells with prescribed density
#ifdef USE_OMP
#pragma omp parallel for schedule(static) num_threads(dom.Nproc)
#endif
for (size_t i=0; i<dat.Right.Size(); ++i)
{
Cell * c = dat.Right[i];
if (c->IsSolid) continue;
double vx = -1.0 + (c->F[0]+c->F[2]+c->F[4] + 2.0*(c->F[1]+c->F[5]+c->F[8]))/c->RhoBC;
c->F[3] = c->F[1] - (2.0/3.0)*c->RhoBC*vx;
c->F[7] = c->F[5] - (1.0/6.0)*c->RhoBC*vx + 0.5*(c->F[2]-c->F[4]);
c->F[6] = c->F[8] - (1.0/6.0)*c->RhoBC*vx - 0.5*(c->F[2]-c->F[4]);
c->Rho = c->VelDen(c->Vel);
}
}
int main(int argc, char **argv) try
{
size_t nproc = 1;
if (argc==2) nproc=atoi(argv[1]);
double u_max = 0.1; // Poiseuille's maximum velocity
double Re = 100000.0; // Reynold's number
size_t nx = 800;
size_t ny = 100;
int radius = ny/10 + 1; // radius of inner circle (obstacle)
double nu = u_max*(2*radius)/Re; // viscocity
LBM::Domain Dom(D2Q9, nu, iVec3_t(nx,ny,1), 1.0, 1.0);
UserData dat;
Dom.UserData = &dat;
dat.vmax = u_max;
//Assigning the left and right cells
for (size_t i=0;i<ny;i++)
{
dat.Left .Push(Dom.Lat[0].GetCell(iVec3_t(0 ,i,0)));
dat.Right.Push(Dom.Lat[0].GetCell(iVec3_t(nx-1,i,0)));
// set parabolic profile
double L = ny - 2; // channel width in cell units
double yp = i - 1.5; // ordinate of cell
double vx = dat.vmax*4/(L*L)*(L*yp - yp*yp); // horizontal velocity
dat.Vel.Push(vx);
dat.Left [i]->RhoBC = 1.0;
//dat.Right[i]->VelBC = 0.08,0.0,0.0;
dat.Right[i]->RhoBC = 1.0;
}
dat.rho = 1.0;
// set inner obstacle
int obsX = ny/2; // x position
int obsY = ny/2+3; // y position
for (size_t i=0;i<nx;i++)
for (size_t j=0;j<ny;j++)
{
if ((i-obsX)*(i-obsX)+(j-obsY)*(j-obsY)<radius*radius)
{
Dom.Lat[0].GetCell(iVec3_t(i,j,0))->IsSolid = true;
}
}
//Assigning solid boundaries at top and bottom
for (size_t i=0;i<nx;i++)
{
Dom.Lat[0].GetCell(iVec3_t(i,0 ,0))->IsSolid = true;
Dom.Lat[0].GetCell(iVec3_t(i,ny-1,0))->IsSolid = true;
}
double rho0 = 1.0;
Vec3_t v0(0.08,0.0,0.0);
//Initializing values
for (size_t i=0;i<Dom.Lat[0].Ncells;i++)
{
Dom.Lat[0].Cells[i]->Initialize(rho0, v0);
}
//Solving
Dom.Time = 0.0;
Dom.Solve(40000.0,80.0,Setup,NULL,"tlbm01",true,nproc);
}
MECHSYS_CATCH
Boundary Conditions
It can be noticed that both an structure UserData and one function
Setup are defined. These two items allow the user to interact with the
internal for loop representing the time evolution inside the function Solve.
The UserData contains the information the user may want to control or report,
in this case it contains five elements. It is worth noting that the boundary conditions of this problem involves a
velocity boundary condition for the fluid entering trough the left limit and a pressure (density) boundary condition
in the right hand size.
struct UserData
{
Array Left;
Array Right;
Array Vel;
double vmax;
double rho;
};
| |
The first array Left will store the pointers to the cells in the left
boundary of the domain, and Right does the same with the right hand size. The
array of doubles Vel will store the velocity profile which in the case of a
pipe flow is the well known parabolic profile. Finally vmax and rho are both the maximum velocity of the profile and the density for the right
hand size boundary.
Boundary conditions will be imposed with the function Setup:
void Setup (LBM::Domain & dom, void * UD)
{
UserData & dat = (*static_cast<UserData *>(UD));
// Cells with prescribed velocity
#ifdef USE_OMP
#pragma omp parallel for schedule(static) num_threads(dom.Nproc)
#endif
for (size_t i=0; i<dat.Left.Size(); ++i)
{
Cell * c = dat.Left[i];
if (c->IsSolid) continue;
double rho = (c->F[0]+c->F[2]+c->F[4] + 2.0*(c->F[3]+c->F[6]+c->F[7]))/(1.0-dat.Vel[i]);
c->F[1] = c->F[3] + (2.0/3.0)*rho*dat.Vel[i];
c->F[5] = c->F[7] + (1.0/6.0)*rho*dat.Vel[i] - 0.5*(c->F[2]-c->F[4]);
c->F[8] = c->F[6] + (1.0/6.0)*rho*dat.Vel[i] + 0.5*(c->F[2]-c->F[4]);
c->Rho = c->VelDen(c->Vel);
}
// Cells with prescribed density
#ifdef USE_OMP
#pragma omp parallel for schedule(static) num_threads(dom.Nproc)
#endif
for (size_t i=0; i<dat.Right.Size(); ++i)
{
Cell * c = dat.Right[i];
if (c->IsSolid) continue;
double vx = -1.0 + (c->F[0]+c->F[2]+c->F[4] + 2.0*(c->F[1]+c->F[5]+c->F[8]))/c->RhoBC;
c->F[3] = c->F[1] - (2.0/3.0)*c->RhoBC*vx;
c->F[7] = c->F[5] - (1.0/6.0)*c->RhoBC*vx + 0.5*(c->F[2]-c->F[4]);
c->F[6] = c->F[8] - (1.0/6.0)*c->RhoBC*vx - 0.5*(c->F[2]-c->F[4]);
c->Rho = c->VelDen(c->Vel);
}
}
The first line of the function will cast the void pointer UD into one of
the clase UserData. The next stage is a for loop going trough all the cells in
the left hand size which is the one with imposed velocity. That for is parallelized by the pragma precompiler command which takes the variable Nproc of
the LBM::Domain class which will be set up later. Inside the loop some
elements are important:
if (c->IsSolid) continue;
double rho = (c->F[0]+c->F[2]+c->F[4] + 2.0*(c->F[3]+c->F[6]+c->F[7]))/(1.0-dat.Vel[i]);
c->F[1] = c->F[3] + (2.0/3.0)*rho*dat.Vel[i];
c->F[5] = c->F[7] + (1.0/6.0)*rho*dat.Vel[i] - 0.5*(c->F[2]-c->F[4]);
c->F[8] = c->F[6] + (1.0/6.0)*rho*dat.Vel[i] + 0.5*(c->F[2]-c->F[4]);
c->Rho = c->VelDen(c->Vel);
The flag IsSolid of the Cell class
is a boolean variable determining if the particle is a solid node or not. If it is solid the boundary condition does
not apply and the bounce-back condition will be imposed. The next lines are the formulas for the distribution
functions F going in the directions 1,5 and 8 which are the only ones with
components going to the right direction. They will ensure that the velocity on those nodes is the same given by the
array dat.Vel which contains the parabolic velocity profile as mentioned
before. These formulas as well as the bounce back condition can be found in the LBM
literature with ease.
Similarly the next for loop goes through the right boundary. The formulas are similar but now they ensure that the
density on those nodes is the ones given by the variable RhoBC.
Cell * c = dat.Right[i];
if (c->IsSolid) continue;
double vx = -1.0 + (c->F[0]+c->F[2]+c->F[4] + 2.0*(c->F[1]+c->F[5]+c->F[8]))/c->RhoBC;
c->F[3] = c->F[1] - (2.0/3.0)*c->RhoBC*vx;
c->F[7] = c->F[5] - (1.0/6.0)*c->RhoBC*vx + 0.5*(c->F[2]-c->F[4]);
c->F[6] = c->F[8] - (1.0/6.0)*c->RhoBC*vx - 0.5*(c->F[2]-c->F[4]);
c->Rho = c->VelDen(c->Vel);
The main scope
The first line of the main scope declares an integer variable nproc which
will be imported from the command line when the program is called.
size_t nproc = 1;
if (argc==2) nproc=atoi(argv[1]);
The next steps assign the maximun velocity of the parabolic velocity profile in the left hand size, the Reynolds
number, the dimensions of the array nx,ny the radius of the obstacle and the
viscosity nu obtained from the Reynolds number and the velocity. Then the LBM
domain is created. The first argument from the constructor is the type of LBM method which in this case is the D2Q9 scheme whose definition can be easily found in the LBM literature. The
second argument is the viscosity of the fluid. The third is a vector of integers containing the lattice dimensions.
The last two arguments are the grid size and time step which are taken as one. Finally an object of the UserData class is declared and assigned to the variable Dom.UserData. This step is necessary for the application of the boundary
conditions explained above.
double u_max = 0.1; // Poiseuille's maximum velocity
double Re = 100000.0; // Reynold's number
size_t nx = 800;
size_t ny = 100;
int radius = ny/10 + 1; // radius of inner circle (obstacle)
double nu = u_max*(2*radius)/Re; // viscocity
LBM::Domain Dom(D2Q9, nu, iVec3_t(nx,ny,1), 1.0, 1.0);
UserData dat;
Dom.UserData = &dat;
Next, the arrays Left, Right and Vel from the UserData. For it the command GetCell is useful to get a pointer to the cell in the coordinates given by the
integer vector serving as argument. Also the density for the cells at the right hand size are imposed.
dat.vmax = u_max;
//Assigning the left and right cells
for (size_t i=0;i<ny;i++)
{
dat.Left .Push(Dom.Lat[0].GetCell(iVec3_t(0 ,i,0)));
dat.Right.Push(Dom.Lat[0].GetCell(iVec3_t(nx-1,i,0)));
// set parabolic profile
double L = ny - 2; // channel width in cell units
double yp = i - 1.5; // ordinate of cell
double vx = dat.vmax*4/(L*L)*(L*yp - yp*yp); // horizontal velocity
dat.Vel.Push(vx);
dat.Left [i]->RhoBC = 1.0;
//dat.Right[i]->VelBC = 0.08,0.0,0.0;
dat.Right[i]->RhoBC = 1.0;
}
dat.rho = 1.0;
The obstacle is obtained by going through the whole domain and setting some cells as solid if they fulfill the circle equation. The obstacle is positioned above the middle line to ensure the appearance of turbulent eddies. Also the upper and bottom boundaries are set as solid as well to model pipe flow.
int obsX = ny/2; // x position
int obsY = ny/2+3; // y position
for (size_t i=0;i<nx;i++)
for (size_t j=0;j<ny;j++)
{
if ((i-obsX)*(i-obsX)+(j-obsY)*(j-obsY)<radius*radius)
{
Dom.Lat[0].GetCell(iVec3_t(i,j,0))->IsSolid = true;
}
}
//Assigning solid boundaries at top and bottom
for (size_t i=0;i<nx;i++)
{
Dom.Lat[0].GetCell(iVec3_t(i,0 ,0))->IsSolid = true;
Dom.Lat[0].GetCell(iVec3_t(i,ny-1,0))->IsSolid = true;
}
Finally an initial velocity and density are assigned for the fluid for all the cells. Then the function Solve is called. This function is very similar as the one called in the example of two polyhedra colliding in the DEM domain. the one difference is that is
does not require a time step since it was declared in the constructor. Again the first argument is the final time and
the second is the timespan for outputs. The third is the name of the function Setup containing the boundary condition formulas. This function is going to be called at each time step.
The fourth argument is a NULL pointer which can be filled with another function
similar to Setup with the only difference being that this one is going to be
called only during each output step. The fifth argument is the file key. The fourth is true if you want to visualize the
results and the last one is the number of processors for parallelazation.
double rho0 = 1.0;
Vec3_t v0(0.08,0.0,0.0);
//Initializing values
for (size_t i=0;i<Dom.Lat[0].Ncells;i++)
{
Dom.Lat[0].Cells[i]->Initialize(rho0, v0);
}
//Solving
Dom.Time = 0.0;
Dom.Solve(40000.0,80.0,Setup,NULL,"tlbm01",true,nproc);
Running the simulation and visualizing the results
After following the compiling instructions the example should be run by typing:
user@mach: ~/msys_tutorials $ ./test_cylinder 8
Where the number 8 is the number of processors. Please check your hardware specifications before running.
Visualizing requires again the Visit software as in the colliding example. The only difference is that now the pseudocolor layer should be chosen for the velocity field. Likewise the velocity as vectors can also be plotted. The velocity field should look like this video: