MercuryDPM  Beta
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Advanced tutorials

Overview

This page contains the following tutorials:

/Drivers/MercurySimpleDemos:

/Drivers/ChuteDemos:

Particles driven by a rotating coil

Problem description

File CoilSelfTest.cpp treats particles placed in a feeder being driven forward by a rotating coil. The entire code of this problem can be viewed here: Particles driven by a rotating coil (code).

Headers

The following headers are included:

These are basically the headers used for the beginner tutorials, except for the additional Coil class.

Class CoilSelfTest

CoilSelfTest, like many of the previous tutorials, inherits from the Mercury3D class.

class CoilSelfTest : public Mercury3D {


The different components of the class will be explained in turn in the following.

Data members

CoilSelfTest has the following (public) data members:

public:
double particleRadius;
Coil* coil;
InfiniteWall* frontWall;
InfiniteWall* backWall;
InfiniteWall* topWall;
InfiniteWall* bottomWall;
InfiniteWall* leftWall;

The six walls are what makes up the rectangular container, which initially contains the particles. The right wall contains a hole, which serves as an exit for the particles. The coil is an object on its own (of the Coil class, obviously), and is to be placed right in front of the container exit hole. particleRadius and species determine the particles' radius and species (i.e., material properties), respectively.

SetupInitialConditions()

The method which sets up the initial conditions of the problem takes up the major part of the driver.
The method starts by a specification of some elementary problem properties, like direction of gravity (-Y), particle radius and the geometrical size of the system.

// gravity, particle radius
setGravity(Vec3D(0.0, -9.8, 0.0));
particleRadius = 0.2;
// set problem geometry
setXMax(1.0);
setYMax(5.0);
setZMax(2.0);
setXMin(-1.0);
setYMin(-1.0);

Next, the particle species is defined. The particles in this problem will use a linear visco-elastic (normal) contact model. The dissipation and stiffness defining the contact model can be set in different ways. In this example the contact model properties are set by giving a collision time, coefficient of restitution and the particle mass.

species = speciesHandler.copyAndAddObject(LinearViscoelasticSpecies());
species->setDensity(1000);
species->setCollisionTimeAndRestitutionCoefficient(0.05, 0.8, pow(particleRadius, 3) * constants::pi * 4.0 / 3.0 * 1000);

The particles are initially contained by a container made up of six walls, five of which are defined to be infinite. The wall in the positive Z-direction is different. Its normal is set (rightWall->set(...) in the positive Z-direction (Vec3D(0, 0, 1)), is set a distance zMax_ (which is returned by getZMax()) from the origin, and contains a 'hole' (which is practically a horizontal tube, since the wall is 'infinite') around the X-axis of radius 1.0 in order to let the particles through.

frontWall = wallHandler.copyAndAddObject(InfiniteWall());
frontWall->set(Vec3D(-1, 0, 0), Vec3D(getXMin(), 0, 0));
backWall = wallHandler.copyAndAddObject(InfiniteWall());
backWall->set(Vec3D(1, 0, 0), Vec3D(getXMax(), 0, 0));
bottomWall = wallHandler.copyAndAddObject(InfiniteWall());
bottomWall->set(Vec3D(0, -1, 0), Vec3D(0, getYMin(), 0));
topWall = wallHandler.copyAndAddObject(InfiniteWall());
topWall->set(Vec3D(0, 1, 0), Vec3D(0, getYMax(), 0));
leftWall = wallHandler.copyAndAddObject(InfiniteWall());
leftWall->set(Vec3D(0, 0, -1), Vec3D(0, 0, getZMin()));
rightWall = wallHandler.copyAndAddObject(InfiniteWallWithHole());
rightWall->set(Vec3D(0, 0, 1), getZMax(), 1.0);

After that the Coil object is added. Its properties are subsequently set by using the Coil::set() method, specifying its starting position (the origin, i.e. Vec3D(0,0,0)), its length (1.0), its radius (1.0 - particleRadius), its number of turns (2.0), its rotation speed (-1.0 revelations per second) and its thickness (0.5 * particleRadius).
NB: its direction or central axis is not specified, since these are the Z-direction and the Z-axis, respectively, by default.

// creation of the coil and setting of its properties
coil = wallHandler.copyAndAddObject(Coil());
// set(Start position, Length, Radius, Number of turns, Rotation speed, Thickness)
coil->set(Vec3D(0, 0, 0), 1.0, 1.0 - particleRadius, 2.0, -1.0, 0.5 * particleRadius);

The particle properties are set subsequently. The particleHandler is cleared just to be sure, the particle to be copied into the container is created and the previously set species is assigned to it. The particle is assigned a zero velocity and the previously defined particle radius.

particleHandler.clear();
p0.setSpecies(species);
p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
p0.setRadius(particleRadius);

After specifying the particle properties, the container is filled with copies of the particle. In this example, particles are placed in a rectangular grid pattern, on evenly spaced positions. First, the number of particles that fit in each direction is computed. Then, a triple for-loop passes every possible particle position, and a particle is placed on the position if there's no part of the coil there.

// Nx*Ny*Nz particles are created and placed on evenly spaced positions in
// the domain [xMin_,xMax_]*[yMin_,yMax_]*[zMin_,zMax_] (unless their position
// is already occupied by the coil).
int Nx = static_cast<int> (std::floor((getXMax() - getXMin()) / (2.1 * particleRadius)));
int Ny = static_cast<int> (std::floor((getYMax() - getYMin()) / (2.1 * particleRadius)));
int Nz = static_cast<int> (std::floor((getZMax() - getZMin()) / (2.1 * particleRadius)));
Mdouble distance;
Vec3D normal;
// Place particles
for (int i = 0; i < Nx; i++)
for (int j = 0; j < Ny; j++)
for (int k = 0; k < Nz; k++) {
p0.setPosition(Vec3D(getXMin()+(getXMax() - getXMin())*(0.5 + i) / Nx, getYMin()+(getYMax() - getYMin())*(0.5 + j) / Ny, getZMin()+(getZMax() - getZMin())*(0.5 + k) / Nz));
if (!coil->getDistanceAndNormal(p0, distance, normal)) {
particleHandler.copyAndAddObject(p0);
} else {
//std::cout<<p0.getPosition()<<std::endl;
}
}

actionsBeforeTimeStep()

The actionsBeforeTimeStep() method specifies all actions that need to be performed at the beginning of each time step, i.e. before the actual numerical solution at that time step is computed. In this case, it states that from 1 second into the simulation time and onward, the coil should turn with the given rotational speed.

void actionsBeforeTimeStep() {
if (getTime() > 1)
coil->move_time(getTimeStep());
}

int main()

In the main program, a CoilSelfTest object is created, after which some of its basic properties are set like its number of dimensions (three), time step and saving parameters. Lastly, the problem is actually solved by calling its solve() method.

int main(int argc UNUSED, char *argv[] UNUSED) {
// create CoilSelfTest object
CoilSelfTest problem;
// set some basic problem properties
problem.setName("CoilSelfTest");
problem.setSystemDimensions(3);
problem.setTimeStep(0.02 * 0.05);
problem.setTimeMax(2.0);
problem.setSaveCount(helpers::getSaveCountFromNumberOfSavesAndTimeMaxAndTimestep(1000, problem.getTimeMax(), problem.getTimeStep()));
// actually solving the problem
problem.solve();
}

(Return to Overview of advanced tutorials)

Particles on an inclined chute

Problem description

File ChuteDemo.cpp treats particles cascading down an inclined chute. The entire code of this problem can be viewed here: Particles on an inclined chute (code).

Headers

The following headers are included:

The particle species is manually set here, and therefore needs to be included. We're treating a chute problem here, and the Chute class needs to be included therefore as well. The particles, walls and boundaries classes are already implemented by the chute class and don't need inclusion here.

int main()

Since the whole structure of the problem is already implemented in the Chute class, no separate class needs to be set up. The setup of initial conditions of the Chute class is shortly treated at the end of this section.
The main driver program starts by initialising a Chute object. Next, the most basic problem properties are set, namely its name (which determines the naming of the data output files), save count (which is the number of time steps skipped between every saved one), particle collision time (which is a species property, but also used in setting the time step), time step and maximum time. Note, that the total number of time steps saved to the output files is not directly set, but is equal to the maximum time divided by time step size and save count.


Next, the particle properties are set. setFixedParticleRadius() sets the radius of the fixed chute bottom particles, while setInflowParticleRadius() sets the inflow particles to be monidisperse with the given particle radius. If inflow particles with random radii are desired, setMinInflowParticleRadius() and setMaxInflowParticleRadius() can be used instead to set the minimum and maximum particle radius, respectively.
The particle species (i.e. its intrinsic material properties) are set next, by specifying the density (species.setDensity()) and the characteristic collision time and coefficient of restitution (with a typical particle mass given; species.setCollisionTimeAndRestitutionCoefficient( \( t_c, r_c, m\))).


The chute properties are subsequently set by specifying the chute's length, width and angle relative to the horizontal.


The chute inflow parameters (besides the previously set inflow particle properties) are set by specifying the inflow height (in Z-direction), the mean iflow particle velocity (in X-direction), and the particle velocity variance (in ratio of the mean velocity).


After all the problem parameters are specified, the simulation is run by calling the solve() method.

problem.solve();

Chute::setupInitialConditions()

{
logger(FATAL, "[Chute::setupInitialConditions()] Chute % cannot complete because no species have been "
"defined.", getName());
// create the chute's side walls in Y-direction (which are solid if the chute is not periodic)
// create a particle of which (altered) copies will fill the chute insertion boundary
p->setSpecies(speciesHandler.getObject(0)); // by default, insert particles of species 0
// set up the insertion boundary and add to handler
//creates the bottom of the chute
}

The setup of initial conditions of the Chute class starts by checking for the presence of a species, and returns an error if there is none. Make sure therefore that you assign a species to the Chute object's speciesHandler before you call the solve() method.
After that, the side walls (in the Y-direction) are set up by calling Chute::setSideWalls(). These are set to be solid, infinite walls by default, but can be set to be periodic instead by setting Chute::isChutePeriodic_ to be true.
A particle is then created (on the heap) which is assigned the first (and only) species in the speciesHandler which we earlier specified in the driver (ChuteDemo.cpp).
A ChuteInsertionBoundary is created, and its parameters subsequently set by its set() method. The set() method arguments are, respectively:

  • p: the previously specified particle
  • maxFailed_: internally used parameter
  • Vec3D(getXMin(), getYMin(), getZMin()): the first defining corner of the cuboidal insertion boundary
  • Vec3D(getXMax(), getYMax(), getZMax()): the second defining corner of the cuboidal insertion boundary
  • min- / maxInflowParticleRadius_: the minimum and maximum radii of inflow particles
  • fixedParticleRadius_: the particle radius making up the chute bottom
  • inflowVelocity(Variance)_: the mean velocity of inflow particles and the allowed variance about the mean

After setting the insertion boundary characteristics, it is added to the problem's boundaryHandler.
Lastly, the chute's bottom is created. The type of bottom created may be set by calling the Chute::setRoughBottomType() method in the driver, giving either of the following four arguments (which are of type enum RoughBottomType):

  • FLAT: just a flat wall (of the species given in the driver)
  • MONOLAYER_ORDERED: (a single layer of) fixed particles in a rectangular grid pattern
  • MONOLAYER_DISORDERED: (default) a single layer of randomly placed particles
  • MULTILAYER: a few layers of randomly placed particles, with a random variation in vertical position as well

(Return to Overview of advanced tutorials)

Particles on a chute with a multilayered bottom

Problem description

File roughBottomDemo.cpp treats particles cascading down an inclined chute with a rough bottom consisting of multiple layers of fixed particles. The entire code of this problem can be viewed here: Particles on a chute with a multilayered bottom (code).

Headers

Class ...

Data members

SetupInitialConditions()

actionsBeforeTimeStep()

int main()

Chute with hopper

Problem description

File roughBottomDemo.cpp treats particles which flow from a hopper onto an inclined chute below. The entire code of this problem can be viewed here: Particles on a chute with a multilayered bottom (code).

Headers

Class ...

Data members

SetupInitialConditions()

actionsBeforeTimeStep()

int main()

Axisymmetric hopper

Problem description

Headers

Class ...

Data members

SetupInitialConditions()

actionsBeforeTimeStep()

int main()