Particles driven by a rotating coil

Problem description:

This tutorial simulates particle flows around a coil for a short time. The code can be found in CoilSelfTest.cpp. It treats particles placed in a feeder being driven forward by a rotating coil.

Headers

The following headers are included:

#include "Mercury3D.h"
#include "Walls/InfiniteWall.h"
#include "Walls/InfiniteWallWithHole.h"
#include "Walls/Coil.h"
#include "Species/LinearViscoelasticSpecies.h"

These are basically the headers used for the beginner tutorials, except for the additional Coil class.

Before the main function

CoilSelfTest, like many of the previous tutorials, inherits from the Mercury3D class.

class CoilSelfTest : public Mercury3D 

The different components of the class will be explained in turn in the following.

step 1: Define Walls
The particles are initially contained by a container made up of six walls, five of which are defined to be infinite. The wall in the positive Z-direction is different. Its normal is set (rightWall->set(...) in the positive Z-direction (Vec3D(0, 0, 1)), is set a distance zMax_ (which is returned by getZMax()) from the origin, and contains a 'hole' (which is practically a horizontal tube, since the wall is 'infinite') around the X-axis of radius 1.0 in order to let the particles through.

        InfiniteWall w;
        w.setSpecies(speciesHandler.getObject(0));
        //front wall
        w.set(Vec3D(-1, 0, 0), Vec3D(getXMin(), 0, 0));
        wallHandler.copyAndAddObject(w);
 
        //back wall
        w.set(Vec3D(1, 0, 0), Vec3D(getXMax(), 0, 0));
        wallHandler.copyAndAddObject(w);
 
        //bottom wall
        w.set(Vec3D(0, -1, 0), Vec3D(0, getYMin(), 0));
        wallHandler.copyAndAddObject(w);
 
        //top wall
        w.set(Vec3D(0, 1, 0), Vec3D(0, getYMax(), 0));
        wallHandler.copyAndAddObject(w);
 
        //left wall
        w.set(Vec3D(0, 0, -1), Vec3D(0, 0, getZMin()));
        wallHandler.copyAndAddObject(w);
 
        //right wall
        InfiniteWallWithHole rightWall;
        rightWall.setSpecies(speciesHandler.getObject(0));
        rightWall.set(Vec3D(0, 0, 1), getZMax(), 1.0);
        wallHandler.copyAndAddObject(rightWall);

step 2: Coil
After that the Coil object is added. Its geometrical properties are subsequently set by using the Coil::set() method, specifying its starting position (the origin, i.e. Vec3D(0,0,0)), its length (1.0), its radius (1.0 - particleRadius), its number of turns (2.0), its rotation speed (-1.0 revelations per second) and its thickness (0.5 * particleRadius).
NB: its direction or central axis is not specified, since these are the Z-direction and the Z-axis, respectively, by default.

        // creation of the coil and setting of its properties
        Coil coil;
        coil.setSpecies(speciesHandler.getObject(0));
        // the syntax to set the coil geometry is as follows:
        // set(Start position, Length, Radius, Number of turns, Rotation speed, Thickness)
        coil.set(Vec3D(0, 0, 0), 1.0, 1.0 - particleRadius, 2.0, -1.0, 0.5 * particleRadius);
        auto pCoil = wallHandler.copyAndAddObject(coil);

step 3: Create Particles
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 previously set species is assigned to it. The particle is assigned a zero velocity and the previously defined particle radius.

        particleHandler.clear();
        SphericalParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
        p0.setRadius(particleRadius);

step 4: Place Particles
After specifying the particle properties, the container is filled with copies of the particle. In this example, particles are placed in a rectangular grid pattern, on evenly spaced positions. First, the number of particles that fit in each direction is computed. Then, a triple for-loop passes every possible particle position, and a particle is placed on the position if there's no part of the coil there.

        // Nx*Ny*Nz particles are created and placed on evenly spaced positions in
        // the domain [xMin_,xMax_]*[yMin_,yMax_]*[zMin_,zMax_] (unless their position
        // is already occupied by the coil).
        const unsigned int Nx = static_cast<unsigned int> (std::floor((getXMax() - getXMin()) / (2.1 * particleRadius)));
        const unsigned int Ny = static_cast<unsigned int> (std::floor((getYMax() - getYMin()) / (2.1 * particleRadius)));
        const unsigned int Nz = static_cast<unsigned int> (std::floor((getZMax() - getZMin()) / (2.1 * particleRadius)));
        Mdouble distance;
        Vec3D normal;
 
        for (unsigned int i = 0; i < Nx; i++)
        {
            for (unsigned int j = 0; j < Ny; j++)
            {
                for (unsigned int k = 0; k < Nz; k++)
                {
                    p0.setPosition(Vec3D(getXMin() + (getXMax() - getXMin()) * (0.5 + i) / Nx,
                                         getYMin() + (getYMax() - getYMin()) * (0.5 + j) / Ny,
                                         getZMin() + (getZMax() - getZMin()) * (0.5 + k) / Nz));
                    if (!pCoil->getDistanceAndNormal(p0, distance, normal)) //if there is no collision with the coil
                    {
                        particleHandler.copyAndAddObject(p0);
                    }
                    else
                    {
                        logger(DEBUG, "particle at position % could not be inserted", p0.getPosition());
                    }
                }
            }
        }

Data members:
CoilSelfTest only has two (public) data members, namely a pointer to the coil:

    Coil* coil;
    Mdouble particleRadius;

Actions Before TimeStep:
The actionsBeforeTimeStep() method specifies all actions that need to be performed at the beginning of each time step, i.e. before the actual numerical solution at that time step is computed. In this case, it states that from 1 second into the simulation time and onward, the coil should turn with the given rotational speed.

    void actionsBeforeTimeStep() override
    {
        if (getTime() > 1)
        {
            coil->move_time(getTimeStep());
        }
    }

Main Function

In the main program, a CoilSelfTest object is created, after which some of its basic properties are set like its number of dimensions (three), time step and saving parameters. Lastly, the problem is actually solved by calling its solve() method.

int main(int argc UNUSED, char* argv[] UNUSED)
{
    
    // create CoilSelfTest object
    CoilSelfTest problem;
 
    // set some basic problem properties
    problem.setName("CoilSelfTest");
    problem.setSystemDimensions(3);
    problem.setGravity(Vec3D(0.0, -9.8, 0.0));
 
    // set problem geometry
    problem.setXMax(1.0);
    problem.setYMax(5.0);
    problem.setZMax(2.0);
    problem.setXMin(-1.0);
    problem.setYMin(-1.0);
    problem.setTimeMax(0.5);
    problem.particleRadius = 0.2;
 
    LinearViscoelasticSpecies species;
    species.setDensity(1000);
    Mdouble tc = 0.05;
    Mdouble restitutionCoefficient = 0.8;
 
    Mdouble particleMass = pow(problem.particleRadius, 3) * constants::pi * 4.0 / 3.0 * species.getDensity();
    species.setCollisionTimeAndRestitutionCoefficient(tc, restitutionCoefficient, particleMass);
    problem.speciesHandler.copyAndAddObject(species);
 
    problem.setTimeStep(0.02 * 0.05);
    problem.setSaveCount(helpers::getSaveCountFromNumberOfSavesAndTimeMaxAndTimeStep(1000, problem.getTimeMax(),
                                                                                     problem.getTimeStep()));
    problem.solve();
    
}

Next, the particle species is defined. The particles in this problem will 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 the contact model properties are set by giving a collision time, coefficient of restitution and the particle mass.

    LinearViscoelasticSpecies species;
    species.setDensity(1000);
    Mdouble tc = 0.05;
    Mdouble restitutionCoefficient = 0.8;
 
    Mdouble particleMass = pow(problem.particleRadius, 3) * constants::pi * 4.0 / 3.0 * species.getDensity();
    species.setCollisionTimeAndRestitutionCoefficient(tc, restitutionCoefficient, particleMass);
    problem.speciesHandler.copyAndAddObject(species);

Reference:

Imole, O. I., Krijgsman, D., Weinhart, T., Magnanimo, V., Montes, B. E. C., Ramaioli, M., & Luding, S. (2016). Reprint of" Experiments and discrete element simulation of the dosing of cohesive powders in a simplified geometry". Powder technology, 293, 69-81.

(Return to Overview of advanced tutorials)