Free cooling granular system (3D)

Model/Simulation:

The homogeneous and inhomogeneous regimes of a free cooling granular system (3D) are simulated. Just like the Free cooling granular system (2D), particles are initially randomized with a large kinetic energy and collisions occur dissipatively. This model is common in the study of the kinetic theory of granular gases [1-3].

Below, a short paraview animation of the 3D free cooling demo is shown.


Simulation set up:

In the animation, the simulation volume consists of a box with equal side length and periodic boundary conditions in three dimensions (3D). Just like the 2D demo case (Free cooling granular system (2D)), the initial state with random particle positions and velocities is prepared in the following way:

  1. The particles first sit on a regular lattice and get a random velocity with a total momentum of zero.
  2. Then the simulation is started without dissipation and runs for a reasonable number of collisions per particle so that the system becomes homogeneous, "and the velocity distribution approaches a Maxwellian."
  3. The above state is then used as the initial configuration for the dissipative regimes.

The FreeCooling3DDemo class inherits from the Mercury3D class.

[FCD_3D:headers]
Definition: FreeCooling3DDemo.cpp:43
This adds on the hierarchical grid code for 3D problems.
Definition: Mercury3D.h:37

Thus, the following headers are included:

Data members of the class:\n

The FreeCooling3DDemo class has, but not limited to two (public) data members, namely a pointer to the species and the number of particles:

The key components of the class are explained in-turn, in the following:

step 1: Define Walls/Periodic Boundaries
For experimental purposes, particles can be contained by a three dimensional box made up of six infinite walls. If this is required, then please see the Free cooling granular system (3D) in walls. Otherwise, the simulation volume is replace by periodic boundaries.

pb.set(Vec3D(1,0,0), getXMin(), getXMax());
boundaryHandler.copyAndAddObject(pb);
pb.set(Vec3D(0,1,0), getYMin(), getYMax());
boundaryHandler.copyAndAddObject(pb);
pb.set(Vec3D(0,0,1), getZMin(), getZMax());
boundaryHandler.copyAndAddObject(pb);
Defines a pair of periodic walls. Inherits from BaseBoundary.
Definition: PeriodicBoundary.h:41
void set(Vec3D normal, Mdouble distanceLeft, Mdouble distanceRight)
Defines a PeriodicBoundary by its normal and positions.
Definition: PeriodicBoundary.cc:84
Definition: Vector.h:51

step 2: Create Particles
Next, the particle species is defined. The particles in this problem 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 these contact model parameters are defined.

species.setDensity(2e3);
species.setDissipation(0.0);
species.setStiffness(1e3);
problem.FC3D_Species = species;
problem.speciesHandler.copyAndAddObject(species);
void setDissipation(Mdouble dissipation)
Allows the normal dissipation to be changed.
Definition: LinearViscoelasticNormalSpecies.cc:117
void setStiffness(Mdouble new_k)
Allows the spring constant to be changed.
Definition: LinearViscoelasticNormalSpecies.cc:93
void setDensity(Mdouble density)
Definition: ParticleSpecies.cc:108

The particle properties are set subsequently. The particleHandler is cleared just to be sure it is empty, then the particle to be copied into the container is created and the set species is assigned to it.

particleHandler.clear();
p0.setSpecies(speciesHandler.getObject(0));
void setSpecies(const ParticleSpecies *species)
Definition: BaseParticle.cc:818
A spherical particle is the most simple particle used in MercuryDPM.
Definition: SphericalParticle.h:37

step 3: Place Particles
After specifying the particle properties, the container is filled with copies of the particle. In this example, particles are placed in a lattice grid pattern, on evenly spaced positions.

for (unsigned int i = 0; i < N; ++i)
{
const unsigned int ix = (i % N1);
const unsigned int iz = static_cast<unsigned int>( static_cast<double>(i) / N1 / N1);
const unsigned int iy = (i - ix - N1 * N1 * iz) / N1;
// set particle position
const double x = (getXMax() - getXMin()) * (ix + 1) / (N1 + 1);
const double y = (getYMax() - getYMin()) * (iy + 1) / (N1 + 1);
const double z = (getZMax() - getZMin()) * (iz + 1) / (N1 + 1);
p0.setPosition(Vec3D(x, y, z));
// set random velocities for the particle
p0.setVelocity(Vec3D(random.getRandomNumber(-2.0,2.0), random.getRandomNumber(-2.0,2.0), random.getRandomNumber(-2.0,2.0)));
p0.setRadius(0.0025);
particleHandler.copyAndAddObject(p0);
}
void setVelocity(const Vec3D &velocity)
set the velocity of the BaseInteractable.
Definition: BaseInteractable.cc:350
virtual void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
Definition: BaseInteractable.h:239
virtual void setRadius(Mdouble radius)
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species)
Definition: BaseParticle.cc:553
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51

step 4: Centre of mass velocity
Next, the center of mass velocity is subtracted to ensure a reduced random velocity. This results into a center of mass velocity nearly equal to zero.

// Compute the center of mass velocity
double particle_mass = p0.getMass();
double M_b = N*particle_mass; // mass of the bulk system
Vec3D V_com = {0,0,0};
for (int k = 0; k < particleHandler.getNumberOfObjects() ; k++){
BaseParticle* p = particleHandler.getObject(k);
V_com += (particle_mass*p->getVelocity())/M_b;
}
// Compute the reduced velocity for each particle
for (int k = 0; k < particleHandler.getNumberOfObjects() ; k++){
BaseParticle* p = particleHandler.getObject(k);
p->setVelocity(p->getVelocity() - V_com);
}
virtual const Vec3D & getVelocity() const
Returns the velocity of this interactable.
Definition: BaseInteractable.cc:329
Definition: BaseParticle.h:54
Mdouble getMass() const
Returns the particle's mass.
Definition: BaseParticle.h:322

Actions After TimeStep:
The actionsAfterTimeStep() method specifies all actions that need to be performed in between time steps, i.e. Since, the simulation started with zero dissipation, after a reasonable number of collisions per particle, the system will be homogeneous, and the actionsAfterTimeStep() is used to set the initial configuration for the "next" dissipative regime.

void actionsAfterTimeStep() override{
if (getTime() > 4e5*getTimeStep()) {
FC3D_Species.setDissipation(0.232);
}
}

Main Function

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

int main(int argc UNUSED, char* argv[] UNUSED)
{
// Problem setup
species.setDensity(2e3);
species.setDissipation(0.0);
species.setStiffness(1e3);
problem.FC3D_Species = species;
problem.speciesHandler.copyAndAddObject(species);
problem.N = 1000;
problem.setName("FreeCooling3DDemo");
problem.setGravity(Vec3D(0.0, 0.0, 0.0));
problem.setTimeStep(5e-5);
problem.setSaveCount(4000);
problem.setTimeMax(50.0);
problem.setMax(0.064,0.064,0.064);
problem.setHGridMaxLevels(1);
problem.setParticlesWriteVTK(true);
problem.solve();
}
@ ONE_FILE
all data will be written into/ read from a single file called name_
#define UNUSED
Definition: GeneralDefine.h:39
int main(int argc, char *argv[])
Definition: T_protectiveWall.cpp:215
std::enable_if<!std::is_pointer< U >::value, U * >::type copyAndAddObject(const U &object)
Creates a copy of a Object and adds it to the BaseHandler.
Definition: BaseHandler.h:379
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_ for all files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:408
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. LinearViscoelasticSpecies, etc.
Definition: DPMBase.h:1427
void setName(const std::string &name)
Allows to set the name of all the files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:422
void setFileType(FileType fileType)
Sets File::fileType_ for all files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:459
void setParticlesWriteVTK(bool writeParticlesVTK)
Sets whether particles are written in a VTK file.
Definition: DPMBase.cc:942
void setTimeStep(Mdouble newDt)
Sets a new value for the simulation time step.
Definition: DPMBase.cc:1234
void setTimeMax(Mdouble newTMax)
Sets a new value for the maximum simulation duration.
Definition: DPMBase.cc:873
void solve()
The work horse of the code.
Definition: DPMBase.cc:4270
void setMax(const Vec3D &max)
Sets the maximum coordinates of the problem domain.
Definition: DPMBase.cc:1082
void setGravity(Vec3D newGravity)
Sets a new value for the gravitational acceleration.
Definition: DPMBase.cc:1383
unsigned int N
[FCD_3D:aftertime]
Definition: FreeCooling3DDemo.cpp:109
LinearViscoelasticSpecies FC3D_Species
Definition: FreeCooling3DDemo.cpp:110
void setHGridUpdateEachTimeStep(bool updateEachTimeStep)
Sets whether or not the HGrid must be updated every time step.
Definition: MercuryBase.cc:176
void setHGridMaxLevels(unsigned int HGridMaxLevels)
Sets the maximum number of levels of the HGrid in this MercuryBase.
Definition: MercuryBase.cc:476
void setHGridCellOverSizeRatio(Mdouble cellOverSizeRatio)
Sets the ratio of the smallest cell over the smallest particle.
Definition: MercuryBase.cc:463

References:

  1. Luding, S. (2005). Structure and cluster formation in granular media. Pramana, 64(6), 893-902.
  2. Brilliantov, N. V., & Pöschel, T. (2010). Kinetic theory of granular gases. Oxford University Press.
  3. Pöschel, T., & Luding, S. (Eds.). (2001). Granular gases (Vol. 564). Springer Science & Business Media

(Return to Overview of advanced tutorials)