Five Particles

Problem description:

This tutorial presents five particles positioned on a irregular base made from five fixed particles. The system is nondimensionalised such that particle diameter d = 1, particle mass m = 1 and gravity g = 1. The code can be found in FiveParticles.cpp.

Snapshot of the final state of the Five Particles simulation (left). Coarse-graining is applied to obtain the bulk density ρ (centre) and pressure p (right).

Headers

Before the main function

class FiveParticles : public Mercury3D
{
public:
void setupInitialConditions() override {
//set parameters to define the species properties
const Mdouble particleRadius=0.5; // such that diameter is one
const Mdouble collisionTime = 0.005; //relatively stiff particles
const Mdouble restitution = 0.88; //restitution close to glass particles
//define species
species->setDensity(6.0/constants::pi); //such that particle mass is one
species->setCollisionTimeAndNormalAndTangentialRestitutionCoefficient(collisionTime, restitution, restitution, species->getMassFromRadius(particleRadius)); //set stiffness and dissipation
species->setSlidingFrictionCoefficient(0.5);
//set gravity, time step
setGravity(Vec3D(0,0,-1)); // such that gravity is one
setTimeStep(0.02 * collisionTime);
//set domain size
setXMin(0); setYMin(0); setZMin(0);
setXMax(5); setYMax(1); setZMax(2.5);
//define common particle properties
p0.setRadius(particleRadius);
//define five fixed particles
p0.setPosition(Vec3D(0.5,0.5,0.1));
p0.setPosition(Vec3D(1.5,0.5,-.2));
p0.setPosition(Vec3D(2.5,0.5,0.0));
p0.setPosition(Vec3D(3.5,0.5,0.1));
p0.setPosition(Vec3D(4.5,0.5,0.05));
//define five free particles
p0.unfix();
p0.setPosition(Vec3D(1.,0.5,1.0));
p0.setPosition(Vec3D(2.0,0.5,1.0));
p0.setPosition(Vec3D(3.0,0.5,1.0));
p0.setPosition(Vec3D(4.0,0.5,1.0));
p0.setPosition(Vec3D(1.5,0.5,2.0));
}
};
double Mdouble
Definition: GeneralDefine.h:34
Species< LinearViscoelasticNormalSpecies, SlidingFrictionSpecies > LinearViscoelasticSlidingFrictionSpecies
Definition: LinearViscoelasticSlidingFrictionSpecies.h:34
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
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:613
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
void fixParticle()
Fix Particle function. It fixes a Particle by setting its inverse mass and inertia and velocities to ...
Definition: BaseParticle.cc:167
void unfix()
Unfix Particle function, which required a reference to the Species vector. It unfixes a Particle by c...
Definition: BaseParticle.cc:322
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
void setSpecies(const ParticleSpecies *species)
Definition: BaseParticle.cc:818
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. LinearViscoelasticSpecies, etc.
Definition: DPMBase.h:1427
void setYMin(Mdouble newYMin)
Sets the value of YMin, the lower bound of the problem domain in the y-direction.
Definition: DPMBase.cc:1034
void setYMax(Mdouble newYMax)
Sets the value of YMax, the upper bound of the problem domain in the y-direction.
Definition: DPMBase.cc:1191
void setZMin(Mdouble newZMin)
Sets the value of ZMin, the lower bound of the problem domain in the z-direction.
Definition: DPMBase.cc:1058
void setXMax(Mdouble newXMax)
Sets the value of XMax, the upper bound of the problem domain in the x-direction.
Definition: DPMBase.cc:1165
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created.
Definition: DPMBase.h:1437
void setZMax(Mdouble newZMax)
Sets the value of ZMax, the upper bound of the problem domain in the z-direction.
Definition: DPMBase.cc:1217
RNG random
This is a random generator, often used for setting up the initial conditions etc.....
Definition: DPMBase.h:1432
void setTimeStep(Mdouble newDt)
Sets a new value for the simulation time step.
Definition: DPMBase.cc:1234
void setXMin(Mdouble newXMin)
Sets the value of XMin, the lower bound of the problem domain in the x-direction.
Definition: DPMBase.cc:1010
void setGravity(Vec3D newGravity)
Sets a new value for the gravitational acceleration.
Definition: DPMBase.cc:1383
[FP:headers]
Definition: FiveParticles.cpp:35
void setupInitialConditions() override
This function allows to set the initial conditions for our problem to be solved, by default particle ...
Definition: FiveParticles.cpp:37
This adds on the hierarchical grid code for 3D problems.
Definition: Mercury3D.h:37
Mdouble getRandomNumber()
This is a random generating routine can be used for initial positions.
Definition: RNG.cc:143
A spherical particle is the most simple particle used in MercuryDPM.
Definition: SphericalParticle.h:37
Definition: Vector.h:51
const Mdouble pi
Definition: ExtendedMath.h:45

Main function

int main(int argc, char *argv[])
{
//instantiate the class
FiveParticles fiveParticles;
//set name
fiveParticles.setName("FiveParticles");
//set output and time stepping properties
fiveParticles.setTimeMax(20); //run until the situation is static
fiveParticles.setSaveCount(std::numeric_limits<unsigned >::max()); //save only the first and last time step
//solve
fiveParticles.solve(argc,argv);
logger(INFO,"Execute 'source FiveParticles.sh' to get coarse-grained statistics of the last time step");
helpers::writeToFile("FiveParticles.sh","../MercuryCG/fstatistics FiveParticles -stattype XZ -w 0.1 -h 0.05 -tmin 1 -tmax 30");
logger(INFO,"Run 'FiveParticles.m' in MATLAB/octave to visualise the statistical output");
helpers::writeToFile("FiveParticles.m","addpath('../../MercuryCG/')\n"
"data = loadstatistics('FiveParticles.stat');\n"
"colormap(1-gray)\n"
"contourf(data.x,data.z,data.Density,20,'EdgeColor','none')\n"
"c = colorbar\n"
"c.Label.String = '\\rho';\n"
"title('Density')\n"
"xlabel('x')\n"
"ylabel('z');\n"
"axis equal\n"
"%%\n"
"particles=importdata('FiveParticles.data',' ',12);\n"
"x=particles.data(:,1);\n"
"z=particles.data(:,3);\n"
"r=particles.data(:,7);\n"
"a=linspace(0,2*pi,40);\n"
"xCircle = sin(a);\n"
"zCircle = cos(a);\n"
"hold on;\n"
"for i=1:length(x)\n"
" plot(x(i)+r(i)*xCircle,z(i)+r(i)*zCircle,'Color',.8*[1 1 1])\n"
"end\n"
"hold off");
}
LL< Log::INFO > INFO
Info log level.
Definition: Logger.cc:55
Logger< MERCURYDPM_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
int main(int argc, char *argv[])
Definition: T_protectiveWall.cpp:215
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_ for all files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:408
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 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
bool writeToFile(std::string filename, std::string filecontent)
Writes a string to a file.
Definition: FileIOHelpers.cc:58

Reference:

Weinhart, T., Thornton, A. R., Luding, S., & Bokhove, O. (2012). From discrete particles to continuum fields near a boundary. Granular Matter, 14(2), 289-294.

(Return to Overview of advanced tutorials)