PeriDEM - High-fidelity modeling of granular media consisting of deformable complex-shaped particles

Table of contents

• Introduction
• Documentation
• Tutorial
• Examples
• Brief implementation details
• Installation
  • Dependencies
  • Building the code
  • Recommendations for quick build
• Running simulations
  • Two-particle with wall
  • Compressive test
  • Attrition tests
  • Single particle deformation
• Visualizing results
• Developers

Introduction

Implementation of the high-fidelity model of granular media that combines the advantages of peridynamics and discrete element method (DEM). The model has the following advantages over existing mechanical models for granular media: - handle intra-particle deformation and breakage/damage - handle arbitrary shape of the particle. Inter-particle contact is not specific to any shape of the particle - tunable inter-particle contact parameters - easy to add different mechanical constitutive laws within peridynamics for individual particle deformation

For more details about the model and results, we refer to the paper:

Prashant K. Jha, Prathamesh S. Desai, Debdeep Bhattacharya, Robert P Lipton (2020). Peridynamics-based discrete element method (PeriDEM) model of granular systems involving breakage of arbitrarily shaped particles. Journal of the Mechanics and Physics of Solids, 151, p.104376. Doi https://doi.org/10.1016/j.jmps.2021.104376.

Download pdf here.

We have created channels on various platforms: - PeriDEM on Gitter * Gitter is absolutely open and easy to join. - PeriDEM on slack * Email us if interested in joining the workspace.

Documentation

Doxygen generated documentation details functions and objects in the library.

Tutorial

We explain the setting-up of simulations in further details in tutorial. We consider two-particle test setup with non-circular particles and compressive-test to discuss the various aspects of simulations.

Examples

We next highlight some key examples. Further details are available in examples/README.md.

Two-particle tests

| | | |:-------------------------------------------------------------------------:|:--------------------------------------------------------------------:| | Circular without damping | Circular with damping |

| | | | |:--------------------------------------------------------------------------:|:---------------------------------------------------------------------------:|:------------------------------------------------------------------------------------------------------:| | Different materials | Different radius | Different radius different material |

Two-particle with wall test

| | |:------------------------------------------------------------------------------------------:| | Concave particles |

Compressive tests

Setup for this test consists of 502 circular and hexagonal-shaped particles of varying radius and orientation inside a rectangle container. The container's top wall is moving downward at a prescribed speed, resulting in the compression of the particle system. The quantity of interest is the compressive strength of the media. The reaction force (downward) on the moving wall should increase with the increasing penetration of this wall; however, after a certain amount of compression of the media, the damage will initiate in individual particles, especially those connected by force chains, resulting in the yielding of the system. For more details, we refer to Jha et al. 2021

| | |:-----------------------------------------------------:| | Compressive test simulation |

Attrition tests

We consider mix of different particles in a rotating container. Particles considered include circular, triangular, hexagonal, and drum shaped. Particles come in large and small shapes (their sizes are purturbed randomly). In order to to introduce diversity of material properties, we considered large particles to be tougher compared to the smaller ones. Setup files are in examples/PeriDEM/attrition_tests

| | | |:----------------------------------------------------------------------------------------------------------------------------------------:|:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------:| | Rotating cylinder (setup) | Rotating cylinder with center of rotation offset (setup) |

Complex container geometries can be considered as well. For example, the image below is from attrition_tests and includes rotating rectangle with opening and internal groves of different shapes. The rotating container with particles inside is contained within another rectangle which is fixed in its place.

Single particle deformation

We can use PeriDEM executable or Peridynamics executable in apps directory to simulate the deformation of single particle/structure using peridynamics. See examples/README.md and examples/Peridynamics folder.

Brief implementation details

The main implementation of the model is carried out in the model directory dem. The model is implemented in class DEMModel. Function DEMModel::run() performs the simulation. We next look at some key methods in DEMModel in more details:

DEMModel::run()

This function does three tasks: ```cpp void model::DEMModel::run(inp::Input *deck) { // initialize data init();

// check for restart
if (d_modelDeck_p->d_isRestartActive)
  restart(deck);

// integrate in time
integrate();

} ```

In DEMModel::init(), the simulation is prepared by reading the input files (such as .yaml, .msh, particle_locations.csv files).

DEMModel::integrate()

Key steps in DEMModel::integrate() are ```cpp void model::DEMModel::run(inp::Input *deck) { // apply initial condition if (d_n == 0) applyInitialCondition();

// apply loading
computeExternalDisplacementBC();
computeForces();

// time step
for (size_t i = d_n; i < d_modelDeck_p->d_Nt; i++) {
  // advance simulation to next step
  integrateStep();

  // perform output if needed
  output();
}

} ```

In DEMModel::integrateStep(), we either utilize the central-difference scheme, implemented in DEMModel::integrateCD(), or the velocity-verlet scheme, implemented in DEMModel::integrateVerlet(). As an example, we look at DEMModel::integrateCD() method below: ```cpp void model::DEMModel::integrateVerlet() { // update current position, displacement, and velocity of nodes { tf::Executor executor(util::parallel::getNThreads()); tf::Taskflow taskflow;

  taskflow.for_each_index(
    (std::size_t) 0, d_fPdCompNodes.size(), (std::size_t) 1,
      [this, dt, dim](std::size_t II) {
        auto i = this->d_fPdCompNodes[II];

        const auto rho = this->getDensity(i);
        const auto &fix = this->d_fix[i];

        for (int dof = 0; dof < dim; dof++) {
          if (util::methods::isFree(fix, dof)) {
            this->d_v[i][dof] += 0.5 * (dt / rho) * this->d_f[i][dof];
            this->d_u[i][dof] += dt * this->d_v[i][dof];
            this->d_x[i][dof] += dt * this->d_v[i][dof];
          }
        }
      } // loop over nodes
  ); // for_each

  executor.run(taskflow).get();
}

// advance time
d_n++;
d_time += dt;

// update displacement bc
computeExternalDisplacementBC();

// compute force
computeForces();

// update velocity of nodes (similar to the above) 

} ```

DEMModel::computeForces()

The key method in time integration is DEMModel::computeForces() In this function, we compute internal and external forces at each node of a particle and also account for the external boundary conditions. This function looks like ```cpp void model::DEMModel::computeForces() { // update the point cloud (make sure that dx is updated along with displacment) auto ptcloudupdatetime = dnsearchp->updatePointCloud(dx, true); ptcloudupdatetime += dnsearchp->setInputCloud();

// reset forces to zero ...

// compute peridynamic forces
computePeridynamicForces();

// compute contact forces between particles
computeContactForces();

// Compute external forces
computeExternalForces();

} ```

Further reading

Above gives the basic idea of simulation steps. For more thorough understanding of the implementation, interested readers can look at demModel.cpp.

Installation

Dependencies

Core dependencies are: - cmake (>= 3.10.2) - vtk (>= 7.1.1) - yaml-cpp (>= 0.5.2) - metis (>= 5.1.0) - MPI

Following dependencies are included in the PeriDEM library in external folder (see external/README.md for more details): - fast-cpp-csv-parser (version included - master) - fmt (>= 7.1.3, version included - 10.2.1) - nanoflann (>= 1.3.2, version included - v1.5.5) - taskflow (>= 3.7.0) - doxygen-awesome-css (>= v2.3.3)

Building the code

If all the dependencies are installed on the global path (e.g., /usr/local/), commands for building the PeriDEM code is as simple as ```sh cmake -DEnableDocumentation=OFF -DEnableTests=ON \ -DEnableHighLoadTests=OFF \ -DDisableDockerMPITests=ON \ -DVTKDIR="${VTKDIR}" \ -DMETISDIR="${METISDIR}" \ -DCMAKEBUILDTYPE=Release \

make -j 4 ```

:exclamation: cmake and make commands should be run inside the build directory. You can create the build directory either inside or outside the repository.

:exclamation: As of now, we can only build the library and not install it. This means you will have to build the package in the build directory, and use the build/bin/PeriDEM executable. We plan to provide the method to install the library in the future.

We refer to tools/README.md for further details about installing dependencies and building the library in different ubuntu releases.

Future plans

We are trying to make PeriDEM MPI-friendly so that we can target large problems. We are moving in following key directions: - MPI parallelism for Peridynamics simulations (deformation of single body subjected to external loading) - MPI parallelism for PeriDEM simulations. Issue is distributing particles to different processors and performing communication efficiently - Asynchronous parallelism within MPI? Currently, we use Taskflow to perform parallel for loops in a non-mpi simulation. In the future, we will be interested in using multithreading combined with MPI to further speed-up the simulations - GPU parallelism?

We are looking for collaborators and HPC experts in making the most of available compute resource and performing truly large-scale high-fidelity granular media simulations. If any of the above future directions interest you or if you have new directions in mind, please do reach out to us.

Ask for help

In the past, PeriDEM library depended on large libraries such as HPX, PCL, Boost (explicitly dependence). We have put a lot of efforts into reducing the dependencies to absolutely minimum so that it is easier to build and run PeriDEM in different operating systems and clusters. At this point, only major library it depends on is VTK which can be compiled to different machines quite successfully (patience is needed in compiling VTK though). If you carefully read the information and use the scripts provided, you should be able to compile PeriDEM in ubuntu (>= 18.04) and mac.

Feel free to reach out or open an issue. For more open discussion of issues and ideas, contact via PeriDEM on Gitter or PeriDEM on slack (for slack, email us to join). If you like some help, want to contribute, extend the code, or discuss new ideas, please do reach out to us.

Running simulations

Assuming that the input file is input.yaml and all other files such as .msh file for particle/wall and particle locations file are created and their filenames with paths are correctly provided in input.yaml, we will run the problem (using 4 threads) sh <path of PeriDEM>/bin/PeriDEM -i input.yaml -nThreads 4

Some examples are listed below.

Two-particle with wall

Navigate to the example directory examples/PeriDEM/twoparticles/twopwallconcavediffmaterialdiff_size/inp and run the example as follows ```sh

manually

cd examples/PeriDEM/twoparticles/twopwallconcavediffmaterialdiffsize/inp mkdir ../out # <-- make directory for simulation output. In .yaml, we specify output path as './out' bin/PeriDEM -i input0.yaml -nThreads 2

or call run.sh script

cd examples/PeriDEM/twoparticles/twopwallconcavediffmaterialdiff_size ./run.sh 4 # with 4 threads ```

You may also use the included problem_setup.py to modify simulation parameters and run the simulation using run.sh.

:exclamation: You may need to modify the path of PeriDEM executable in run.sh file.

In all problem_setup.py files in the example and test directory, the main function is create_input_file(). Here we set all model parameters, create .yaml input file, and .geo files for meshing.

Compressive test

Navigate to the example directory examples/PeriDEM/compressivetest/compressionlarge_set/inp and run the example as follows (note that this is a computationally expensive example) ```sh cd examples/PeriDEM/compressivetest/compressionlargeset/inp mkdir ../out bin/PeriDEM -i input0.yaml -nThreads 12

or you can use the run.sh script in the path examples/PeriDEM/compressivetest/compressionlarge_set/

```

As before: - you can modify problem_setup.py, see create_input_file() method, to change the simulation settings - run the simulation using run.sh.

Visualizing results

Simulation files output_*.vtu can be loaded in either ParaView or VisIt.

By default, in all tests and examples, we only output the particle mesh, i.e., a pair of nodal coordinate and nodal volume, and not the finite element mesh (it can be enabled by setting Perform_FE_Out: true within Output block in the input yaml file). After loading the file in ParaView, the first thing to do is to change the plot type from Surface to Point Gaussian. Next, a couple of things to do are: - Adjust the radius of circle/sphere at the nodes by going to the Properties tab on the left side and change the value of Gaussian Radius - You may also want to choose the field to display. For starter, you could select the Damage_Z variable, a ratio of maximum bond strain in the neighborhood of a node and critical bond strain. When the Damage_Z value is below one at a given node, the deformation in the vicinity of that node is elastic, whereas when the value is above 1, it indicates there is at least one node in the neighborhood which has bond strain above critical strain (meaning the bond between these two nodes is broken) - You may also need to rescale the plot by clicking on the Zoom to Data button in ParaView - Lastly, when the Damage_Z is very high at few nodes, you may want to rescale the data to the range, say [0,2] or [0,10], so that it is easier to identify regions with elastic deformation and region with fracture.

Contributing

We welcome contributions to the code. In Future plans section above, some potential directions are listed. Please fork this repository, make changes, and make a pull request to the source branch.

Citations

If this library was useful in your work, we recommend citing the following article:

Jha, P.K., Desai, P.S., Bhattacharya, D. and Lipton, R., 2021. Peridynamics-based discrete element method (PeriDEM) model of granular systems involving breakage of arbitrarily shaped particles. Journal of the Mechanics and Physics of Solids, 151, p.104376.

You can also cite the PeriDEM using zenodo doi:

Prashant K., J. (2024). Peridynamics-based discrete element method (PeriDEM) model of granular systems. Zenodo. https://doi.org/10.5281/zenodo.13888588

Developers

• Prashant K. Jha (pjha.sci@gmail.com, prashant.jha@sdsmt.edu)