Parameter study using a protective wall

Problem description

Particles are released from a specific height, roll through a slope and then a protective wall retains them. The user controls the number of particles inserted, the height of release, slope angle, grain size. The output data contain the .fstat file, which includes the information of numbers of grains contained by the wall and overflowing the wall, total force on the wall, volume fraction, stress inside the retained material. An illustrative description is presented in the figure below.

Schematic representation of the parameter study including a protective wall, h: height of release, l: length, w: width, s: slope angle.

Quick explanation:

There are two ways to run the tutorial. First, input parameters set directly in the main file: ProtectiveWall.cpp. Second, input parameters are read from the command line after the compilation, for instance: ./protectiveWall -Np 500 -r 0.01 -h 0.1 -w 0.25 -l 1.0 -s 15.0 -t 20.0, where -Np: Number of particles, -r: particle radius.

Some initial setups:

Three different examples changing the initial setup.

Simulation:

Particles are released from a specific height, roll through a slope, and they are retained by a protective wall.

Headers

The following headers are included:

#include <Mercury3D.h>
#include "Species/LinearViscoelasticSpecies.h"
#include "Walls/InfiniteWall.h"
#include "Walls/IntersectionOfWalls.h"
#include "Boundaries/CubeInsertionBoundary.h"
#include "Boundaries/CubeDeletionBoundary.h"

These are the headers needed for this tutorial. Further explanation can be found in Beginner tutorials.

Main class:

The main class contains the constructor and member functions.

class protectiveWall : public Mercury3D
{
public:
 
    protectiveWall(int Nump, Mdouble pRadius, Mdouble height, Mdouble width,Mdouble Length, Mdouble slopeAngle)
    {
        setNumParticles = Nump; //Number of particles
        setGlobalRadius = pRadius; //Particle radius
        setParticleHeight = height; //Height of release
        setSlopeAngle = constants::pi / 180.0 * slopeAngle; //Slope angle
 
        wallHandler.setWriteVTK(FileType::ONE_FILE); //For visualization
        setParticlesWriteVTK(true); //For visualization
        setGravity(Vec3D(0.0, 0.0, -9.81)); //Set gravity
 
        Mdouble Hipo = Length/cos(setSlopeAngle); //Slope length
 
        setXMax(Length); //Boundary length
        setYMax(width); //Boundary width
        setZMax(Hipo*sin(setSlopeAngle)); //Boundary height using slope length
 
        slopeSpecies = speciesHandler.copyAndAddObject(LinearViscoelasticSpecies());
        slopeSpecies->setDensity(2500.0); //set the species density
        slopeSpecies->setStiffness(20058.5);//set the spring stiffness.
        slopeSpecies->setDissipation(0.01); //set the dissipation.
 
        particleSpecies = speciesHandler.copyAndAddObject(LinearViscoelasticSpecies());
        particleSpecies->setDensity(2500.0); //set the species density
        particleSpecies->setStiffness(258.5);//set the spring stiffness.
        particleSpecies->setDissipation(0.5); //set the dissipation.
 
        speciesHandler.getMixedObject(0,1)->mixAll(slopeSpecies,slopeSpecies); //Particle-wall interactions
 
        insb = new CubeInsertionBoundary;
        SphericalParticle particles;
        particles.setSpecies(particleSpecies);
    
        insb->set(
                &particles,
                1,
                Vec3D(0.9 * getXMax(), 0.45 * getYMax(), getZMax() + 0.97 * setParticleHeight), //Minimum position
                Vec3D(getXMax(), 0.55 * getYMax(), getZMax() + setParticleHeight), //Maximum position
                Vec3D(0, 0, 0), //Minimum velocity
                Vec3D(0, 0, 0));
        insb = boundaryHandler.copyAndAddObject(insb);
 
        //To delete particles
        delb = new DeletionBoundary;
        delb->set(Vec3D(-1,0,0), getXMin());
        delb =  boundaryHandler.copyAndAddObject(delb);
    }
 
 
    void setupInitialConditions() override
    {
        InfiniteWall slope,lateralwall1, lateralwall2, lateralwall3;
 
        slope.setSpecies(slopeSpecies);
        slope.setPosition(Vec3D(getXMin(),getYMin(),getZMin()));
        slope.setOrientation(Vec3D(sin(setSlopeAngle),0,-cos(setSlopeAngle)));
        wallHandler.copyAndAddObject(slope);
 
        lateralwall1.setSpecies(slopeSpecies);
        lateralwall1.set(Vec3D(0.0, 1.0, 0.0), Vec3D(0.0,getYMax(),0.0));
        wallHandler.copyAndAddObject(lateralwall1);
 
        lateralwall2.setSpecies(slopeSpecies);
        lateralwall2.set(Vec3D(0.0, -1.0, 0.0), Vec3D(0.0,getYMin(),0.0));
        wallHandler.copyAndAddObject(lateralwall2);
 
        IntersectionOfWalls roarWall, proctWall;
        //Rear wall.
        roarWall.addObject(Vec3D(1.0, 0.0, 0.0), Vec3D(getXMax(), 0.0, 0.0));
        roarWall.addObject(Vec3D(0.0, 0.0, -1.0), Vec3D(getXMax(), 0.0, getZMax()+setParticleHeight));
        wallHandler.copyAndAddObject(roarWall);
 
 
        //Protective wall
        heightProtWall = setGlobalRadius*10.0;
        proctWall.setSpecies(slopeSpecies);
        proctWall.addObject(Vec3D(-1.0, 0.0, 0.0), Vec3D(getXMin(), 0.0, 0.0));
        proctWall.addObject(Vec3D(0.0, 0.0, -1.0), Vec3D(getXMin(), 0.0, heightProtWall));
        wallHandler.copyAndAddObject(proctWall);
 
    }
 
    void actionsAfterTimeStep() override {
 
        vol_inserted = insb->getVolumeOfParticlesInserted();
        partDel = delb->getNumberOfParticlesDeleted();
 
        //Specific number of particles
        setVolume = 4.0/3.0 * constants::pi * (setGlobalRadius*setGlobalRadius*setGlobalRadius);
        setVolume *= setNumParticles;
 
        //Stop inserting particles
        if(vol_inserted  >= setVolume && !removed_insb){
            boundaryHandler.removeObject(insb->getIndex());
            removed_insb = true;
            setVolume = vol_inserted;
        }
 
        //Wall force and pressure
        wallForce = fabs(wallHandler.getObject(4)->getForce().X);
        wallPressure = wallForce / (getYMax()*heightProtWall);
    }
 
    //Criterion to stop the simulation, otherwise the simulation stops at maxTime.
    bool continueSolve() const override
    {
        static unsigned int counter = 0;
        if (++counter>100)
        {
            counter=0;
            if (getKineticEnergy()<0.001*getElasticEnergy())
                return false;
        }
        return true;
    }
 
    //Print varaibles in the console/terminal
    void printTime()const override
    {
        logger(INFO, "t=%3, "
                     "tmax=%3, "
                     "# Particles inserted=%3, "
                     "# Particles deleted=%3, "
                     "Volume inserted=%3.6, "
                     "WallForce=%3.6, "
                     "WallPressure=%3.6",
               getTime(), getTimeMax(), setNumParticles - partDel, partDel, vol_inserted, wallForce,
               wallPressure);
    }
 
    //Write variables in the fstat file
    void writeFstatHeader(std::ostream &os) const override {
        os<< getTime() //Current time
          << " " << getTimeMax() //MaxTime
          << " " << setNumParticles - partDel //Particles inserted
          << " " << partDel //Partilcles deleted
          << " " << vol_inserted //Volume inserted
          << " " << wallForce //Wall force
          << " " <<  wallPressure //Wall pressure
          << std::endl;
    }
 
 
private:
    Mdouble setGlobalRadius = 1.0; //By default
    Mdouble setParticleHeight = 0.1; //By default
    Mdouble setSlopeAngle = 20.0; //By default
    Mdouble heightProtWall = 1.0; //By default
 
    CubeInsertionBoundary* insb;
    DeletionBoundary* delb;
    int partDel = 0;
    int setNumParticles;
    Mdouble wallForce = 0.0;
    Mdouble wallPressure = 0.0;
    Mdouble vol_inserted = 0.0;
    Mdouble setVolume = 0.0;
    bool removed_insb = false;
 
    LinearViscoelasticSpecies* slopeSpecies;
    LinearViscoelasticSpecies* particleSpecies;
};
 

Main function

int main(int argc, char* argv[])
{
    //Helper
    logger(INFO,
           "Write in the terminal after the compilation'./protectiveWall -Np 500 -r 0.01 -h 0.1 -w 0.25 -l 1.0 -s 15.0 -t "
           "20.0' to run the program");
    
    int Nump = helpers::readFromCommandLine(argc, argv, "-Np", 50); //50 particles
    Mdouble pRadius = helpers::readFromCommandLine(argc, argv, "-r", 0.01); //0.01 [m]
    Mdouble height = helpers::readFromCommandLine(argc, argv, "-h", 0.1); //0.1 [m]
    Mdouble width = helpers::readFromCommandLine(argc, argv, "-w", 0.25);  //0.25 [m]
    Mdouble length = helpers::readFromCommandLine(argc, argv, "-l", 1.0); //1.0 [m]
    Mdouble slopeAngle = helpers::readFromCommandLine(argc, argv, "-s", 15.0); //15 grades
    
    protectiveWall problem(Nump, pRadius, height, width, length, slopeAngle); //Object
    
    Mdouble simTime = helpers::readFromCommandLine(argc, argv, "-t", 5.0); // 5.0 [s]
    problem.setTimeMax(simTime);
    std::string simName = helpers::readFromCommandLine(argc,argv,"-name",std::string("protectiveWall"));
    problem.setName(simName);
    
    problem.setSaveCount(400);
    problem.setTimeStep(0.005 / 50.0); // (collision time)/50.0
    problem.removeOldFiles();
    problem.solve();
    return 0;
}