Free cooling granular system (3D) in walls

Description:

This tutorial is used to simulate homogeneous and inhomogeneous regimes of a free cooling granular system/particles contained in a cubic box (six infinite walls, InfiniteWall). For detailed descriptions about the free cooling granular systems see references [1-3].

Here, 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 in 3D. Just like other free cooling demo tutorials, 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 FreeCooling3DinWallsDemo class inherits from the Mercury3D class.

class FreeCooling3DinWallsDemo : public Mercury3D 

Thus, the following headers are included:

#include <iostream>
#include <Species/LinearViscoelasticSpecies.h>
#include <Walls/InfiniteWall.h>
#include "Mercury3D.h"
#include "MercuryTime.h"

Data members of the class:\n

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

    unsigned int N;
    LinearViscoelasticSpecies FC3D_Species;

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

step 1: Define Walls
Particles can be contained by a three dimensional box made up of six infinite walls.

        wallHandler.clear();
        InfiniteWall w0;
        w0.setSpecies(speciesHandler.getObject(0));
        w0.set(Vec3D(-1.0, 0.0, 0.0), Vec3D(getXMin(), 0, 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(1.0, 0.0, 0.0), Vec3D(getXMax(), 0, 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(0.0, -1.0, 0.0), Vec3D(0, getYMin(), 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(0.0, 1.0, 0.0), Vec3D(0, getYMax(), 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(0.0, 0.0, -1.0), Vec3D(0, 0, getZMin()));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(0.0, 0.0, 1.0), Vec3D(0, 0, getZMax()));
        wallHandler.copyAndAddObject(w0);

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.

    LinearViscoelasticSpecies species;
    species.setDensity(2e3);
    species.setDissipation(0.0);
    species.setStiffness(1e3);
    problem.FC3D_Species = species;
    problem.speciesHandler.copyAndAddObject(species);

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();
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));

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);
        }

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);
        }

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 FreeCooling3DinWallsDemo 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
    FreeCooling3DinWallsDemo problem;
    LinearViscoelasticSpecies species;
    species.setDensity(2e3);
    species.setDissipation(0.0);
    species.setStiffness(1e3);
    problem.FC3D_Species = species;
    problem.speciesHandler.copyAndAddObject(species);
    problem.N = 1000;
    problem.setName("FreeCooling3DinWallsDemo");
    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.setHGridCellOverSizeRatio(1.2);
    problem.setHGridUpdateEachTimeStep(false);
    problem.setFileType(FileType::ONE_FILE);
    problem.setParticlesWriteVTK(true);
    problem.solve();
}

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)