Simple shear (constant stress) of a cuboidal REV

Problem description

Like the previous tutorials, this one demonstrates how the StressStrainControlBoundary works on the simple shear deformation mode. Here in the example we keep a constant shear rate in x direction while try to maintain a constant stress in y direction. We generate a loose initial packing by randomly inserting the particles into the cuboid, then tell the boundary that we need to keep constant shear rate in x direction with a velocity gradient in y direction. We also tell the boundary that we need to have stress in y direction to reach 100 (the actual unit depends on your scaling of other parameters, if you use SI-unit, i.e. second, metre, kilogram, for all the other parameters, then the stress will be Pascal) and be kept there at close as possible to the target value during the shearing. Note that in the current settings, the particles are following affine movements in a strainrate control. You could also set the strainrate control to false, then the boundary will drag the particles to move.

Here is 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, they are the same as the other REV tutorials.

Before the main function

Like the previous example of isotropic compression, simple shear also 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.05) / 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.05;   //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 = 1.0*1.0*1.0 = 1.0, if failed by 100 times, it stops
        insertion.set(&particle, 100, getMin(), getMax(), Vec3D(0, 0, 0), Vec3D(0, 0, 0));
        insertion.setInitialVolume(1.0);
        insertion.checkBoundaryBeforeTimeStep(this);
        logger(INFO,"Inserted % particles",particleHandler.getSize());
 
        //Set up the deformation mode using the boundaries that are based on user's input
        boundaryHandler.clear();                  // Delete all exist boundaries
        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 (this is associated with the stress control and could be tunned accordingly), 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.

Note that for tuning the gainFactor, if the value is too large, then the stress control correction magnitude is too large, resulting in over correction, or if it is too small, then the correction per timestep is too small then you never reach the target stress value. This gain factor is also associated with the shear rate, because the faster you shear, the higher stress difference is created per timestep, then the corresponding correction on the stress has also to be large. However, one has to bear in mind that if the shear deformation is in similar magnitude as the stress control magnitude, i.e. you move the particle in both shear (x) and stress control (y) directions in similar displacement per timestep, it might cause resonance effect resulting in the oscillation in the stress field.

In the current example, we have set the isStrainRateControlled = true, this means we move particles everytime step using affine movement according to the strainrate tensor. This is because we want to shear everything in a more homogenous way and it will reach steady state faster. If you choose to use boundary movement dragging the particles, then simply set this boolean to false.

The current setup also moves the boundary with its strainrate control, therefore the volume is not conserving. If user would like to have a constant volume shearing, simply set the target stress tensor to zero (this case set target stress in y direction to 0).

//Here we define all the control parameters and solve the problem
int main(int argc UNUSED, char* argv[] UNUSED)
{
    //Here we want to set stress in y direction to constant, 100 is just an example,
    //if you would like to keep the volume constant, then everything here should be set to zero
    Matrix3D stressGoal;
    stressGoal.XX = 0.0;
    stressGoal.YY = 100.0;
    stressGoal.ZZ = 0.0;
    stressGoal.XY = 0.0;
 
    //Here we set the simple shear strainrate tensor, negative sign means compression
    Matrix3D strainRate;
    strainRate.XX = 0.0;
    strainRate.YY = 0.0;
    strainRate.ZZ = 0.0;
    strainRate.XY = 1.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_SimpleShear");
    //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.
Shi, H., Roy, S., Weinhart, T., Magnanimo, V., & Luding, S. (2020). Steady state rheology of homogeneous and inhomogeneous cohesive granular materials. Granular matter, 22(1), 14.

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)