Particles driven by a rotating coil

Problem description

File CoilSelfTest.cpp treats particles placed in a feeder being driven forward by a rotating coil. The entire code of this problem can be viewed here: Particles driven by a rotating coil (code).

Headers

The following headers are included:

#include "Mercury3D.h"
#include "Particles/BaseParticle.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.

Class CoilSelfTest

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.

Data members

CoilSelfTest only has one (public) data member, namely a pointer to the coil:

    Coil* coil;

SetupInitialConditions()

The method which sets up the initial conditions of the problem takes up the major part of the driver.
The method starts by a specification of some elementary problem properties, like direction of gravity (-Y), particle radius and the geometrical size of the system.

        // gravity, particle radius
        setGravity(Vec3D(0.0, -9.8, 0.0));
        const double particleRadius = 0.2;
        
        // set problem geometry
        setXMax(1.0);
        setYMax(5.0);
        setZMax(2.0);
        setXMin(-1.0);
        setYMin(-1.0);

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);
        const double tc = 0.05;
        const double restitutionCoefficient = 0.8;
        const double particleMass = pow(particleRadius, 3) * constants::pi * 4.0 / 3.0 * species.getDensity();
        species.setCollisionTimeAndRestitutionCoefficient(tc, restitutionCoefficient, particleMass);
        speciesHandler.copyAndAddObject(species);

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);

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 = wallHandler.copyAndAddObject(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);

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();
        BaseParticle p0;
        p0.setSpecies(speciesHandler.getObject(0));
        p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
        p0.setRadius(particleRadius);

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;
        
        // Place particles
        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 (!coil->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());
                    }
                }
            }
        }

actionsBeforeTimeStep()

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());
        }

    }

int main()

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.setTimeStep(0.02 * 0.05);
    problem.setTimeMax(0.5);
    problem.setSaveCount(helpers::getSaveCountFromNumberOfSavesAndTimeMaxAndTimestep(1000, problem.getTimeMax(),
                                                                                     problem.getTimeStep()));

    // actually solving the problem
    problem.solve();
    
}

(Return to Overview of advanced tutorials)