DEM Examples
Collision of 2 particles
This tutorial aims to build the simulation for 2 bodies colliding in an elastic way. The user declares a sphero-cube and a sphero-tetrahedron in a collision course. The energy, liner and angular momenta can be measured to check these conservation laws. The whole simulation is described in the following code which will be explained in detail:
// MechSys
#include <mechsys/dem/domain.h>
int main(int argc, char **argv) try
{
DEM::Domain dom;
// add cube
Vec3_t x(-10,0,0);
Vec3_t v(1.,0,0);
Vec3_t w(0,M_PI/5,0);
dom.AddCube (/*Tag*/-1, /*Position*/x,/*SpheroRadius*/0.3,/*Length*/3.,/*Density*/1.);
dom.Particles[0]->v = v;
dom.Particles[0]->w = w;
// add tetrahedron
x = 10.0, 0 , 0;
w = 0.0, 0 , M_PI/10.;
v = -1.0, 0 , 0;
dom.AddTetra (/*Tag*/-1, /*Position*/x,/*SpheroRadius*/0.5,/*Length*/5.,/*Density*/1.);
dom.Particles[1]->v = v;
dom.Particles[1]->w = w;
// particle parameters
Dict D;
D.Set(-1,"Gn Gt Mu",0.0,0.0,0.0);
dom.SetProps(D);
// initial constants
Vec3_t l0;
Vec3_t p0;
double Ek0,Ep0,E0;
dom.LinearMomentum (p0);
dom.AngularMomentum (l0);
E0 = d.CalcEnergy (Ek0,Ep0);
// solve
dom.CamPos = 0.0,30.0,0.0;
dom.Solve(/*tf*/30.0, /*dt*/1.0e-4, /*dtOut*/0.3, NULL, NULL, /*filekey*/"test_dynamics");
// final constants
Vec3_t l1; // initial linear momentum
Vec3_t p1; // initial angular momentum
double Ek1,Ep1,E1; // initial energy
dom.LinearMomentum (p1);
dom.AngularMomentum (l1);
E1 = dom.CalcEnergy (Ek1,Ep1);
// check
double tol = 1.0e-3;
double err_l = norm(l1-l0);
double err_p = norm(p1-p0);
double err_E = fabs(E1-E0);
double error = err_l + err_p + err_E;
cout << "Error in energy = " << err_E << endl;
cout << "Error in angular momentum = " << err_l << endl;
cout << "Error in linear momentum = " << err_p << endl;
// results
if (error > tol) return 1;
else return 0;
}
MECHSYS_CATCH
Initialization
MechSys/DEM simulations are prepared with a C++ source code. The first step is to include the required MechSys libraries. In this case the file domain.h contains the DEM domain definitions:
#include <mechsys/dem/domain.h>
Simulation script structure
The simulation script structure is:
#include ...
int main(int argc, char **argv) try
{
...
}
MECHSYS_CATCH
After including the required MechSys' libraries (and C++ standard libraries, if necessary), the main
function is declared.
For MechSys, it is important to add the try keyword just after
the int main() and before the first opening brace {. Then, after the last closing
brace the MECHSYS_CATCH is added.
Declaring the DEM Domain and Particles
The line,
DEM::Domain d;
declares the container dom of the class
DEM::Domain. This class defines the DEM universe and also has many functions
to include particles, some of which will be used in this tutorial.
To declare the cube, some 3D vectors are needed. The Vec3_t data type is used
for fast vectorial operations:
Vec3_t x(-10,0,0);
Vec3_t v(1.,0,0);
Vec3_t w(0,M_PI/5,0);
The vectors x v w
define respectively the cube center position, its linear velocity and its angular velocity. The declared angular velocity ensures that the cube will
spin around the y-axis counter clockwise. Once these three vectors are declared, the cube is introduced into the domain:
dom.AddCube (/*Tag*/-1, /*Position*/x,/*SpheroRadius*/0.3,/*Length*/3.,/*Density*/1.);
This function creates a cube into the memory. The first argument is a Tag which will be useful later on. The second is the position vector declared above. Next is the Sphero-Radius which is the radius of the rounded corners, and is required for the collision law. Ideally this spheroradius should be large enough so the overlapping between the colliding particles is small compared to the particles dimensions. The 4th argument is the length of the cube's edges and finally the density of the cube is also declared.
The method, however, initializes the particle at rest and without rotation. Therefore, to add the vector velocities, internal variables must be accessed:
dom.Particles[0]->v = v;
dom.Particles[0]->w = w;
The object dom.Particles is an array that contains the memory pointers to all the particles created in the domain dom. To access
the first particle, which is the cube, the index 0 must be introduced. Therefore dom.Particles[0] is a pointer to the cube. The internal vector variables v and w are the cube's linear and angular velocity vectors.
In a similar way the tetrahedron is included at the right side of the cube moving towards it:
x = 10.0, 0 , 0;
w = 0.0, 0 , M_PI/10.;
v = -1.0, 0 , 0;
dom.AddTetra (/*Tag*/-1, /*Position*/x,/*SpheroRadius*/0.5,/*Length*/5.,/*Density*/1.);
dom.Particles[1]->v = v;
dom.Particles[1]->w = w;
With both particles properly defined, some collision parameters must be set. MechSys uses the spheropolyhedra method to model the collision of
solid objects. The collision law is defined by two spring constants Kn and Kt; two dissipative constants Gn and Gt; and a friction coefficient μ. Since this simulation deals with an elastic collision, the
two dissipative constants and the friction coefficient must be set to zero. In MechSys this is achieved by the use of dictionaries as follows:
Dict D;
D.Set(-1,"Gn Gt Mu",0.0,0.0,0.0);
dom.SetProps(D);
Where the dictionary D contains the user information on the new parameters. The first argument of the
function D.Set is the particle's tag. Hence different groups of particles can be associated with different
collision parameters. The elastic parameters Kn and Kt (as Kn and Kt in the dictionary string) can also be included into the dictionary, however for this example we will work
with their default values.
Starting the simulation
To check the conservation laws some values must be calculated before the simulation starts:
Vec3_t l0;
Vec3_t p0;
double Ek0,Ep0,E0;
dom.LinearMomentum (p0);
dom.AngularMomentum (l0);
E0 = d.CalcEnergy (Ek0,Ep0);
The vectors l0 and p0 will store the linear and angular momenta given by the functions dom.LinearMomentum and dom.AngularMomentum. The three numbers Ek0, Ep0 and E0 will be respectively the total kinetic energy, the
total potential energy and the total energy as given by the function dom.CalcEnergy. This variable will be kept in memory to be compared with the situation after the simulation.
With the DEM universe defined, the simulation will begin. First, for future visualization, is useful to define the position of a camera:
dom.CamPos = 0.0,30.0,0.0;
Then the system evolves for as long as defined by the user by means of the dom.Solve:
dom.Solve(/*tf*/30.0, /*dt*/1.0e-4, /*dtOut*/0.3, NULL, NULL, /*filekey*/"test_dynamics");
The system will evolve in this case up to tf = 30.0 seconds of simulation time. The next
argument is the integration time step dt, the smaller the more accurate the simulation is but also the more
time consuming. The 3rd argument is the report time step, in other words, the time span between movie frames. The next two arguments are
pointers to functions which are not explained in this simple example. This functions may be called from inside the Solve function, and hence offer
some control over the evolution of the system from outside. Finally, the last argument is a file key which is used to identify the report and
animation files form this simulation.
At the end of the simulation, the same variables that were measured before it are measured again:
Vec3_t l1;
Vec3_t p1;
double Ek1,Ep1,E1;
dom.LinearMomentum (p1);
dom.AngularMomentum (l1);
E1 = d.CalcEnergy (Ek1,Ep1);
And compared with their initial values:
double tol = 1.0e-3;
double err_l = norm(l1-l0);
double err_p = norm(p1-p0);
double err_E = fabs(E1-E0);
double error = err_l + err_p + err_E;
cout << "Error in energy = " << err_E << endl;
cout << "Error in angular momentum = " << err_l << endl;
cout << "Error in linear momentum = " << err_p << endl;
After compiling and running, the simulation should give this results
Error in energy = 0.128839
Error in angular momentum = 0.000122599
Error in linear momentum = 8.30147e-10
Unfortunately the conservation laws are not exactly met in the DEM and therefore there is a difference between the initial and final values. As mentioned before this error can be tuned down with the integration time step
Simulation visualization
MechSys uses VisIt to visualize simulation results. A full tutorial can be found here. After following the instructions the dynamics of the system should resemble the following video: