DEM+LBM Examples
Sphere in fluid and the measurement of Stokes Law
This tutorial aims combines most of the concepts introduced in the previous version with the novel coupling between DEM and LBM. The code will introduce a LBM grid where the fluid is accelerated with a constant body force and a sphere is immersed in it. The force over the sphere can be calculated to obtain the drag coefficient.
struct UserData
{
std::ofstream oss_ss;
Vec3_t acc;
double nu;
double R;
Array<Cell *> xmin;
Array<Cell *> xmax;
};
void Setup(LBM::Domain & dom, void * UD)
{
UserData & dat = (*static_cast<UserData *>(UD));
#ifdef USE_OMP
#pragma omp parallel for schedule(static) num_threads(dom.Nproc)
#endif
for (size_t i=0;i<dom.Lat[0].Ncells;i++)
{
Cell * c = dom.Lat[0].Cells[i];
c->BForcef = c->Rho*dat.acc;
}
}
void Report(LBM::Domain & dom, void * UD)
{
UserData & dat = (*static_cast<UserData *>(UD));
if (dom.idx_out==0)
{
String fs;
fs.Printf("%s_force.res",dom.FileKey.CStr());
dat.oss_ss.open(fs.CStr());
dat.oss_ss << Util::_10_6 << "Time" << Util::_8s << "Fx" << Util::_8s << "Fy" << Util::_8s << "Fz" << Util::_8s << "Tx" << Util::_8s << "Ty" << Util::_8s << "Tz" << Util::_8s << "Vx" << Util::_8s << "Fx" << Util::_8s << "F" << Util::_8s << "Rho" << Util::_8s << "Re" << Util::_8s << "CD" << Util::_8s << "CDsphere \n";
}
if (!dom.Finished)
{
double M = 0.0;
double Vx = 0.0;
size_t nc = 0;
Vec3_t Flux = OrthoSys::O;
for (size_t i=0;i<dom.Lat[0].Ncells;i++)
{
Cell * c = dom.Lat[0].Cells[i];
if (c->IsSolid||c->Gamma>1.0e-8) continue;
Vx += (1.0 - c->Gamma)*c->Vel(0);
Flux += c->Rho*c->Vel;
M += c->Rho;
nc++;
}
Vx /=dom.Lat[0].Ncells;
Flux/=M;
M /=nc;
double CD = 2.0*dom.Particles[0]->F(0)/(Flux(0)*Flux(0)/M*M_PI*dat.R*dat.R);
double Re = 2.0*Flux(0)*dat.R/(M*dat.nu);
double CDt = 24.0/Re + 6.0/(1.0+sqrt(Re)) + 0.4;
if (Flux(0)<1.0e-12)
{
Flux = OrthoSys::O;
CD = 0.0;
Re = 0.0;
}
dat.oss_ss << Util::_10_6 << dom.Time << Util::_8s << dom.Particles[0]->F(0) << Util::_8s << dom.Particles[0]->F(1) << Util::_8s << dom.Particles[0]->F(2) << Util::_8s << dom.Particles[0]->T(0) << Util::_8s << dom.Particles[0]->T(1) << Util::_8s << dom.Particles[0]->T(2) << Util::_8s << Vx << Util::_8s << Flux(0) << Util::_8s << norm(Flux) << Util::_8s << M << Util::_8s << Re << Util::_8s << CD << Util::_8s << CDt << std::endl;
}
else
{
dat.oss_ss.close();
}
}
int main(int argc, char **argv) try
{
if (argc<2) throw new Fatal("This program must be called with one argument: the name of the data input file without the '.inp' suffix.\nExample:\t %s filekey\n",argv[0]);
size_t Nproc = 1;
if (argc==3) Nproc=atoi(argv[2]);
String filekey (argv[1]);
String filename (filekey+".inp");
if (!Util::FileExists(filename)) throw new Fatal("File <%s> not found",filename.CStr());
ifstream infile(filename.CStr());
String ptype;
bool Render = true;
size_t nx = 100;
size_t ny = 50;
size_t nz = 50;
double nu = 0.01;
double dx = 1.0;
double dt = 1.0;
double Dp = 0.1;
double R = 10.0;
double w = 0.001;
double ang= 0.0;
double Tf = 40000.0;
{
infile >> ptype; infile.ignore(200,'\n');
infile >> Render; infile.ignore(200,'\n');
infile >> nx; infile.ignore(200,'\n');
infile >> ny; infile.ignore(200,'\n');
infile >> nz; infile.ignore(200,'\n');
infile >> nu; infile.ignore(200,'\n');
infile >> dx; infile.ignore(200,'\n');
infile >> dt; infile.ignore(200,'\n');
infile >> Dp; infile.ignore(200,'\n');
infile >> R; infile.ignore(200,'\n');
infile >> w; infile.ignore(200,'\n');
infile >> ang; infile.ignore(200,'\n');
infile >> Tf; infile.ignore(200,'\n');
}
LBM::Domain Dom(D3Q15, nu, iVec3_t(nx,ny,nz), dx, dt);
Dom.Step = 1;
Dom.Sc = 0.0;
UserData dat;
Dom.UserData = &dat;
dat.acc = Vec3_t(Dp,0.0,0.0);
dat.R = R;
dat.nu = nu;
for (int i=0;i<nx;i++)
for (int j=0;j<ny;j++)
for (int k=0;k<nz;k++)
{
Vec3_t v0(0.0,0.0,0.0);
Dom.Lat[0].GetCell(iVec3_t(i,j,k))->Initialize(1.0,v0);
}
if (ptype=="sphere") Dom.AddSphere(-1,Vec3_t(0.5*nx*dx,0.5*ny*dx,0.5*nz*dx),R,3.0);
else if (ptype=="tetra" )
{
double e = pow(sqrt(2)*M_PI,1.0/3.0)*2*R;
Dom.AddTetra(-1,Vec3_t(0.5*nx*dx,0.5*ny*dx,0.5*nz*dx),0.05*e,e,3.0,M_PI/4.0,&OrthoSys::e2);
Quaternion_t q;
NormalizeRotation(35.26*M_PI/180.0,OrthoSys::e1,q);
Dom.Particles[0]->Rotate(q,Dom.Particles[0]->x);
NormalizeRotation(ang*M_PI/180.0,OrthoSys::e2,q);
Dom.Particles[0]->Rotate(q,Dom.Particles[0]->x);
}
else if (ptype=="cube" )
{
double e = pow(M_PI/6.0,1.0/3.0)*2*R;
Dom.AddCube(-1,Vec3_t(0.5*nx*dx,0.5*ny*dx,0.5*nz*dx),0.05*e,e,3.0,ang*M_PI/180.0,&OrthoSys::e2);
}
else
{
DEM::Particle * Pa = new DEM::Particle(-1, ptype.CStr(), 0.05*R, 3.0, R );
Dom.Particles.Push(Pa);
Vec3_t t(0.5*nx,0.5*ny,0.5*nz);
Pa->Position(t);
}
Dom.Particles[0]->FixVeloc();
Dom.Particles[0]->w = Vec3_t(0.0,0.0,w);
//Solving
Dom.Solve(Tf,0.01*Tf,Setup,Report,filekey.CStr(),Render,Nproc);
}
MECHSYS_CATCH
UserData Structure
The first step is to define the UserData structure as in the case of the flow past a cylinder
struct UserData
{
std::ofstream oss_ss;
Vec3_t acc;
double nu;
double R;
};
The first variable is an object of the class std::ofstream which is a class
from C++ for output files. The vector acc will contain the acceleration which
is chosen by the user. Part of the output is to report the Reynolds number which requires both the equivalent radius
of the particle R and the viscosity of the fluid nu.
Setup Function
As in the example of flow around cylinder a user defined function is used to internally modified some variables within the time for loop. In this case the function will only applied the body force for each LBM cell providing a constant acceleration.
void Setup(LBM::Domain & dom, void * UD)
{
UserData & dat = (*static_cast<UserData *>(UD));
#ifdef USE_OMP
#pragma omp parallel for schedule(static) num_threads(dom.Nproc)
#endif
for (size_t i=0;i<dom.Lat[0].Ncells;i++)
{
Cell * c = dom.Lat[0].Cells[i];
c->BForcef = c->Rho*dat.acc;
}
}
As can be seen, the for loop modifies the variable BForcef of the Cell class by assigning it a value of the cell density (c->Rho) times the imposed acceleration dat.acc. The #pragma preamble is used for parallelization of
this for loop.
Report Function
Since the goal of this test is to measure the drag coefficient, it is important to report the force over the particle for each time step. The Report function help us to achieve this goal.
void Report(LBM::Domain & dom, void * UD)
{
UserData & dat = (*static_cast<UserData *>(UD));
if (dom.idx_out==0)
{
String fs;
fs.Printf("%s_force.res",dom.FileKey.CStr());
dat.oss_ss.open(fs.CStr());
dat.oss_ss << Util::_10_6 << "Time" << Util::_8s << "Fx" << Util::_8s << "Fy" << Util::_8s << "Fz" << Util::_8s << "Tx" << Util::_8s << "Ty" << Util::_8s << "Tz" << Util::_8s << "Vx" << Util::_8s << "Fx" << Util::_8s << "F" << Util::_8s << "Rho" << Util::_8s << "Re" << Util::_8s << "CD" << Util::_8s << "CDsphere \n";
}
if (!dom.Finished)
{
double M = 0.0;
double Vx = 0.0;
size_t nc = 0;
Vec3_t Flux = OrthoSys::O;
for (size_t i=0;i<dom.Lat[0].Ncells;i++)
{
Cell * c = dom.Lat[0].Cells[i];
if (c->IsSolid||c->Gamma>1.0e-8) continue;
Vx += (1.0 - c->Gamma)*c->Vel(0);
Flux += c->Rho*c->Vel;
M += c->Rho;
nc++;
}
Vx /=dom.Lat[0].Ncells;
Flux/=M;
M /=nc;
double CD = 2.0*dom.Particles[0]->F(0)/(Flux(0)*Flux(0)/M*M_PI*dat.R*dat.R);
double Re = 2.0*Flux(0)*dat.R/(M*dat.nu);
double CDt = 24.0/Re + 6.0/(1.0+sqrt(Re)) + 0.4;
if (Flux(0)<1.0e-12)
{
Flux = OrthoSys::O;
CD = 0.0;
Re = 0.0;
}
dat.oss_ss << Util::_10_6 << dom.Time << Util::_8s << dom.Particles[0]->F(0) << Util::_8s << dom.Particles[0]->F(1) << Util::_8s << dom.Particles[0]->F(2) << Util::_8s << dom.Particles[0]->T(0) << Util::_8s << dom.Particles[0]->T(1) << Util::_8s << dom.Particles[0]->T(2) << Util::_8s << Vx << Util::_8s << Flux(0) << Util::_8s << norm(Flux) << Util::_8s << M << Util::_8s << Re << Util::_8s << CD << Util::_8s << CDt << std::endl;
}
else
{
dat.oss_ss.close();
}
}
The first exception (dom.idx_out==0) checks if the simulation is running
the first instance of the time for loop. It will open the file "%s_force.res" (where %s will be replace by the filekey
introduced in the function Solve). Once opened, the function will write a
header on the file with the tags for the columns storing the time, the components of the force over the particle (F)
, the torque (T), the Reynolds number (Re) and the drag coefficient (CD).
The second exception (!dom.Finished) checks if the last instance of the
time for loop has not been reached (i.e. that if the simulation is still running). The it will go through each cell to
calculate the average velocity. With this, the drag coefficient and the Reynolds number are calculated and stored in
the output file. In particular, the drag coefficient requires to force over the particle which is obtained from the
variable dom.Particles[0]->F. Finally once the last instance of the time for loop is reached, the output file is closed and flushed
out.
The main loop
As in the previous examples, the main loop is used to set the scene for the simulation and call the Solve function. One difference between the previous examples is that now an
input file must be provided which allows to make changes on the parameters without having to recompile to code again.
the input file (test_input.inp) includes all the key variables for the study
sphere = ptype Particle type
0 = Render Produce visualization files?
240 = nx x dimension of lbm grid
60 = ny y dimension of lbm grid
60 = nz z dimension of lbm grid
1.0 = nu viscosity of the fluid
0.4 = dx grid size
1.6e-2 = dt time step
3.1e-2 = Dp fluid acceleration
3.6 = R particle size
0.00 = w angular velocity
45.0 = ang angle of the body around the z axis
1.0e3 = Tf time
Each of this variables are read within the main loop in the following section
infile >> ptype; infile.ignore(200,'\n');
infile >> Render; infile.ignore(200,'\n');
infile >> nx; infile.ignore(200,'\n');
infile >> ny; infile.ignore(200,'\n');
infile >> nz; infile.ignore(200,'\n');
infile >> nu; infile.ignore(200,'\n');
infile >> dx; infile.ignore(200,'\n');
infile >> dt; infile.ignore(200,'\n');
infile >> Dp; infile.ignore(200,'\n');
infile >> R; infile.ignore(200,'\n');
infile >> w; infile.ignore(200,'\n');
infile >> ang; infile.ignore(200,'\n');
infile >> Tf; infile.ignore(200,'\n');
And as such, the input values become variables of the main loop that can be used to define the problem parameters.
For instance, one key variable is the string ptype which determines the
particle type and has four possible options (sphere, tetra, cube and the name of a JSON file to import a mesh
particle). Once the particle is defined, it is added to the domain. For example, for the case of a sphere we have
if (ptype=="sphere") Dom.AddSphere(-1,Vec3_t(0.5*nx*dx,0.5*ny*dx,0.5*nz*dx),R,3.0);
Adds a sphere at the center of the domain with radius given by the variable R and density equal to 3.0. Most of the functions to add particles come from the DEM domain and the
spheropolyhedra method is embedded into the LBM domain as well so that complex shapes particles can interact with each
other and with the fluid.
The last line before the Solve function fix to zero the velocity of the particle
by means of the function FixVeloc and assign an angular velocity along the z
axis. When the particle has its velocities fixed, neither the forces or torques will change the velocity or angular
velocity although they can still be recorded. That ensures that the particle will remain at the center and will keep
spinning.
Dom.Particles[0]->FixVeloc();
Dom.Particles[0]->w = Vec3_t(0.0,0.0,w);
Results
To run the simulation the command ./test_spherefluid test_input 8 must be
entered where 8 can be replaced for any other number of cores. This example is particularly time consumming so it is
advised to use all the available cores. By plotting the results of the file test_input_force.res, a comparison can be
done with empirical formulae found from
experiments for the drag coefficient of a sphere. In the next images, the fluid velocity field and the results for the
lift coefficient for several angular velocities are presented by comparison. The results are taken from this
paper which was written while Mechsys was developed. Please
refer to it for further details.
