HourGlass2DDemo.cpp File Reference

#include "Mercury3D.h"
#include "Particles/BaseParticle.h"
#include "Walls/IntersectionOfWalls.h"
#include "Walls/AxisymmetricIntersectionOfWalls.h"
#include <iostream>
#include "Species/LinearViscoelasticFrictionSpecies.h"

Go to the source code of this file.

Classes
class	HourGlass2D

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

Function Documentation

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

Definition at line 106 of file HourGlass2DDemo.cpp.

References HourGlass2D::ContractionHeight, HourGlass2D::ContractionWidth, BaseHandler< T >::copyAndAddObject(), mathsFunc::cubic(), ParticleSpecies::getDensity(), HourGlass2D::MaxParticleRadius, HourGlass2D::MinParticleRadius, HourGlass2D::N, constants::pi, ParticleSpecies::setDensity(), DPMBase::setGravity(), Files::setName(), Files::setSaveCount(), DPMBase::setTimeMax(), DPMBase::setTimeStep(), DPMBase::setXBallsAdditionalArguments(), DPMBase::setXMax(), DPMBase::setXMin(), DPMBase::setYMax(), DPMBase::setYMin(), DPMBase::setZMax(), DPMBase::setZMin(), DPMBase::solve(), and DPMBase::speciesHandler.

{
    std::cout<< "Hourglass Simulation" <<std::endl;
    // 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());
    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);
    //Calculates collision time for two copies of a particle of given dissipation_, k, effective mass
    std::cout << "MinParticleMass =" << MinParticleMass << std::endl;
    //Calculates collision time for two copies of a particle of given dissipation_, k, effective mass
    Mdouble tc = species->getCollisionTime(MinParticleMass);
    std::cout << "tc  =" << tc << std::endl;
    //Calculates restitution coefficient for two copies of given dissipation_, k, effective mass
    Mdouble r = species->getRestitutionCoefficient(MinParticleMass);
    std::cout << "r   =" << r << std::endl;
    //Calculates the maximum relative velocity allowed for a normal collision of two particles of radius r and particle mass m (for higher velocities particles could pass through each other)
    //std::cout << "vmax=" << helpers::getMaximumVelocity(species->getStiffness(), HGgetSpecies(0)->getDissipation(), HG.MinParticleRadius, MinParticleMass) << std::endl;

    //set other simulation parameters
    HG.setTimeStep(tc / 50.0);
    HG.setTimeMax(3.0);
    HG.setSaveCount(500);
    HG.setXBallsAdditionalArguments("-v0 -solidf");
    HG.N=100; //number of particles
    std::cout << "N   =" << HG.N << std::endl;

    HG.solve(argc, argv);
    return 0;
}