Isotropic compression of a cuboidal REV

Problem description

This tutorial presents how to use the StressStrainControlBoundary to achieve isotropic compression in a 3D cuboidal representative element volume. The particle bed is initially inserted into the box in a loose packing and then isotropically compressed using a constant strainrate control.

Here a short video of the code's output.

Headers

The following headers are included:

#include <Mercury3D.h>
#include <Species/LinearViscoelasticSpecies.h>
#include <Boundaries/StressStrainControlBoundary.h>
#include <Boundaries/CubeInsertionBoundary.h>

These are the necessary headers to be included.

Before the main function

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

class StressStrainControl : public Mercury3D 

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

Step 1: Define the constructor:
Here we first define the constructor that contains user inputs (target stress, strainrate, gainFactor and the switch Boolean key for strainrate control on particle movements) which will be inserted later on in the main function. Then we define here also the domain size and specie.

public:
    //Create the class and passing the user inputs using the constructors
    StressStrainControl(const Matrix3D& stressGoal, const Matrix3D& strainRate, const Matrix3D& gainFactor,
                        bool isStrainRateControlled)
            : stressGoal_(stressGoal), strainRate_(strainRate), gainFactor_(gainFactor),
              isStrainRateControlled_(isStrainRateControlled)
    {
        //Define the domain size
        setMin(Vec3D(0, 0, 0));
        setMax(Vec3D(1, 1, 1));
        //Calculate the mass from smallest particle
        double mass = 2000 * constants::pi * mathsFunc::cubic(2.0 * 0.08) / 6.0;
        //Define the species for the particles
        LinearViscoelasticSpecies species;
        species.setDensity(2000);
        species.setStiffnessAndRestitutionCoefficient(10000, 0.4, mass);
        speciesHandler.copyAndAddObject(species);
    }

Step 2: Set up initial conditions:
Here we setup the system dimensions, timestepping and gravity, then add particles using the CubeInsertionBoundary. After we inserted 96 particles, we then setup the StressStrainControlBoundary using user inputs. For the current case, the final boundary consists of 3 Periodic boundaries in x, y, z directions, respectively.

private:
    //Here we set up the system parameters, adding particles and define boundaries
    void setupInitialConditions() override
    {
        //Set up the micro parameters.
        double rhop = 2000;     //particle density
        double K1 = 10000;      //particle stiffness
        double en = 0.4;        //restitution coefficient
        double radius = 0.08;   //particle radius
        //Calculate the mass from smallest particle
        double mass = rhop * constants::pi * mathsFunc::cubic(2.0 * radius) / 6.0;
        //Calculate the contact duration
        double tc = std::sqrt(mass / 2.0 / K1 * (mathsFunc::square(constants::pi) +
                                                 mathsFunc::square(log(en))));
        setSystemDimensions(3); //set the 3d simulation
        setTimeStep(tc / 50); //set the timestep based on the shortest contact duration
        setGravity(Vec3D(0.0, 0.0, 0.0)); //set gravity in the system
 
        //Add particles
        CubeInsertionBoundary insertion;
        SphericalParticle particle(speciesHandler.getObject(0));
        //Insert 96 particles in the subvolume = 0.8*0.8*0.8 = 0.512, if failed by 100 times, it stops
        insertion.set(&particle, 100, 0.8 * getMin(), 0.8 * getMax(), Vec3D(0, 0, 0), Vec3D(0, 0, 0));
        insertion.setInitialVolume(1.0);
        insertion.checkBoundaryBeforeTimeStep(this);
        logger(INFO,"Inserted % particles",particleHandler.getSize());
 
        //Delete all existing boundaries
        boundaryHandler.clear();
        //Set up the deformation mode using the boundaries that are based on user's input
        StressStrainControlBoundary boundary;
        boundary.setHandler(&boundaryHandler);
        boundary.set(stressGoal_, strainRate_, gainFactor_, isStrainRateControlled_);
        boundaryHandler.copyAndAddObject(boundary);
 
        
    }

Step 3: Declare private data members:
REV_IsotropicCompression has three Matrix3D variables and one Boolean variable. These variables are used to pass the user inputs via the constructor.

    //initialize the private variables for passing the user's inputs in the constructor
    Matrix3D stressGoal_;
    Matrix3D strainRate_;
    Matrix3D gainFactor_;
    bool isStrainRateControlled_;

Main Function

In the main function, we first define the target stress tensor, strainrate tensor, gainFactor tensor (only used if stress tensor is non-zero), and the boolean variable isStrainRateControlled.
Then a StressStrainControl object is created, after which some of its basic properties are set like its name, saveCount and outputs ON/OFF. Finally, the problem is actually solved by calling its solve() method.

//Here we define all the control parameters and solve the problem
int main(int argc UNUSED, char* argv[] UNUSED)
{
    //We want strainrate control, so here we set the target Stress to free, all to zero
    Matrix3D stressGoal;
    stressGoal.XX = 0.0;
    stressGoal.YY = 0.0;
    stressGoal.ZZ = 0.0;
    stressGoal.XY = 0.0;
 
    //Here we set the isotropic strainrate tensor, negative sign means compression
    Matrix3D strainRate;
    strainRate.XX = -0.02;
    strainRate.YY = -0.02;
    strainRate.ZZ = -0.02;
    strainRate.XY = 0.0;
 
    //This is default gainFactor, if you choose to use stress control, you might want to tune this one
    Matrix3D gainFactor;
    gainFactor.XY = 0.0001;
    gainFactor.XX = 0.0001;
    gainFactor.YY = 0.0001;
    gainFactor.ZZ = 0.0001;
 
    //This is by default set to true, so particles are controlled by affine movement.
    //If set it to false, then the particles will only follow the boundary movement.
    bool isStrainRateControlled = true;
    
    //Define the object and problem to solve
    StressStrainControl dpm(stressGoal, strainRate, gainFactor, isStrainRateControlled);
    
    //Set name
    dpm.setName("Tutorial_IsotropicCompression");
    //Set output and time stepping properties
    dpm.setTimeMax(10);
    //Set the SaveCount, i.e. 100 timesteps saving one snapshot
    dpm.setSaveCount(100);
    //Currently all the normal file outputs are switched off, simply switch it on by replacing NO_FILE to ONE_FILE
//    dpm.dataFile.setFileType(FileType::NO_FILE);
//    dpm.restartFile.setFileType(FileType::NO_FILE);
//    dpm.fStatFile.setFileType(FileType::NO_FILE);
//    dpm.eneFile.setFileType(FileType::NO_FILE);
    //Set output the particle information in VTK for ParaView Visualisation
    dpm.setParticlesWriteVTK(true);
    //Because of periodic boundary, out put wall files is not necessary in this case
    //dpm.wallHandler.setWriteVTK(true);
    //Solve the problem
    dpm.solve();
}

Reference:

Imole, O. I., Kumar, N., Magnanimo, V., & Luding, S. (2013). Hydrostatic and shear behavior of frictionless granular assemblies under different deformation conditions. KONA Powder and Particle Journal, 30, 84-108.

If you would like to know more about the definitions used in the StressStrainControlBoundary and various deformation possibilities, please click here (right click and open in a new tab if left click does not work): StressStrainControlBoundary.pdf

(Return to Overview of advanced tutorials)