HourGlass2DDemo.cpp File Reference

#include <iostream>
#include "Mercury3D.h"
#include "Walls/IntersectionOfWalls.h"
#include "Species/LinearViscoelasticFrictionSpecies.h"
#include <cstdlib>
#include <chrono>

Classes
class	HourGlass2D

Functions
int	main (int argc, char *argv[])

Function Documentation

main()

int main	(	int argc	,
		char *	argv[]
	)

{
    // Start measuring elapsed time
    std::chrono::time_point<std::chrono::system_clock> startClock, endClock;
    startClock = std::chrono::system_clock::now();
    
    logger(INFO, "Hourglass Simulation");
    // note: this code is based on stefan's implementation, see
    // /storage2/usr/people/sluding/COMPUTERS/hpc01/sluding/MDCC/MDCLR/DEMO/W7
    // however, it is scaled to SI units by the scaling factors
    // d=1e-3, m=1e-6, g=1
    
    //all parameters should be set in the main file
    //here, we use SI units (generally, other, scaled units are possible)
    
    //create an instance of the class and name it
    HourGlass2D HG;
    HG.setName("HourGlass2D");
    LinearViscoelasticFrictionSpecies* species = HG.speciesHandler.copyAndAddObject(LinearViscoelasticFrictionSpecies());
    //LinearViscoelasticSpecies* species = HG.speciesHandler.copyAndAddObject(LinearViscoelasticSpecies());
    species->setDensity(2000);
 
    //specify geometry
    //specify dimensions of the hourglass
    Mdouble Width = 10e-2; // 10cm
    Mdouble Height = 60e-2; // 60cm
    //specify how big the wedge of the contraction should be
    Mdouble ContractionWidth = 2.5e-2; //2.5cm
    Mdouble ContractionHeight = 5e-2; //5cm
    //set domain accordingly (domain boundaries are not walls!)
    HG.setXMin(0.0);
    HG.setXMax(Width);
    HG.setYMin(0.0);
    HG.setYMax(Width);
    HG.setZMin(0.0);
    HG.setZMax(Height);
    //these parameters are needed in setupInitialConditions()
    HG.ContractionWidth = ContractionWidth;
    HG.ContractionHeight = ContractionHeight;
 
    //specify particle properties
    species->setDensity(2000.0);
    //these parameters are needed in setupInitialConditions()
    HG.MinParticleRadius = 6e-3; // 6mm
    HG.MaxParticleRadius = 10e-3; //10mm
 
    //specify body forces
    HG.setGravity(Vec3D(0.0, 0.0, -9.8));
 
    //specify contact properties
    //normal forces
    species->setStiffness(1e5);
    species->setDissipation(0.63);
    //tangential (sliding) forces
    species->setSlidingFrictionCoefficient(0.5);
    species->setSlidingStiffness(1.2e4);
    species->setSlidingDissipation(0.16);
    //tangential (rolling) torques
    species->setRollingFrictionCoefficient(0.2);
    species->setRollingStiffness(1.2e4);
    species->setRollingDissipation(6.3e-2);
    //normal (torsion/spin) torques
    species->setTorsionFrictionCoefficient(0.1);
    species->setTorsionStiffness(1.2e4);
    species->setSlidingDissipation(6.3e-2);
    
    //test normal forces
    Mdouble MinParticleMass =
            species->getDensity() * 4.0 / 3.0 * constants::pi * mathsFunc::cubic(HG.MinParticleRadius);
    logger(INFO, "MinParticleMass =%\n", MinParticleMass, Flusher::NO_FLUSH);
    Mdouble tc = species->getCollisionTime(MinParticleMass);
    logger(INFO, "tc  =%\n"
                 "r   =%\n"
                 "vmax=%\n",
           tc, species->getRestitutionCoefficient(MinParticleMass),
           species->getMaximumVelocity(HG.MinParticleRadius, MinParticleMass), Flusher::NO_FLUSH);
    
    //set other simulation parameters
    HG.setTimeStep(tc / 50.0);
    HG.setTimeMax(2.0);//run until 3.0 to see full simulation
    HG.setSaveCount(500);
    HG.setXBallsAdditionalArguments("-v0 -solidf");
    //uncomment next two line 2 to create paraview files
    HG.setParticlesWriteVTK(true);
    HG.wallHandler.setWriteVTK(FileType::ONE_FILE);
    
    HG.N = 100; //number of particles
    logger(INFO, "N   = %", HG.N);
 
    HG.solve(argc, argv);
    
    // Measure elapsed time
    endClock = std::chrono::system_clock::now();
    std::chrono::duration<double> elapsed_seconds = endClock - startClock;
    logger(INFO, "Elapsed time for solving the PDE: % s", elapsed_seconds.count());
    
    return 0;
}

References HourGlass2D::ContractionHeight, HourGlass2D::ContractionWidth, BaseHandler< T >::copyAndAddObject(), mathsFunc::cubic(), e(), LinearViscoelasticNormalSpecies::getCollisionTime(), ParticleSpecies::getDensity(), LinearViscoelasticNormalSpecies::getMaximumVelocity(), LinearViscoelasticNormalSpecies::getRestitutionCoefficient(), GlobalParameters::Height, INFO, logger, HourGlass2D::MaxParticleRadius, HourGlass2D::MinParticleRadius, HourGlass2D::N, NO_FLUSH, ONE_FILE, constants::pi, ParticleSpecies::setDensity(), LinearViscoelasticNormalSpecies::setDissipation(), DPMBase::setGravity(), DPMBase::setName(), DPMBase::setParticlesWriteVTK(), FrictionSpecies::setRollingDissipation(), FrictionSpecies::setRollingFrictionCoefficient(), FrictionSpecies::setRollingStiffness(), DPMBase::setSaveCount(), SlidingFrictionSpecies::setSlidingDissipation(), SlidingFrictionSpecies::setSlidingFrictionCoefficient(), SlidingFrictionSpecies::setSlidingStiffness(), LinearViscoelasticNormalSpecies::setStiffness(), DPMBase::setTimeMax(), DPMBase::setTimeStep(), FrictionSpecies::setTorsionFrictionCoefficient(), FrictionSpecies::setTorsionStiffness(), WallHandler::setWriteVTK(), DPMBase::setXBallsAdditionalArguments(), DPMBase::setXMax(), DPMBase::setXMin(), DPMBase::setYMax(), DPMBase::setYMin(), DPMBase::setZMax(), DPMBase::setZMin(), DPMBase::solve(), DPMBase::speciesHandler, and DPMBase::wallHandler.