Modelling large deformations with
Material Point Method


Krishna Kumar, kks32@cam.ac.uk
University of Cambridge.

Kenichi Soga
University of California, Berkeley.




Cambridge-Berkeley Computational Geomechanics

  • Lattice-Boltzmann + Discrete Element Method
  • Finite Element Method - Thermo-Hydro Mechanical Coupling
  • Material Point Method
  • Lattice Element Method
View the CB-Geo website for more information and software tools

LBM-DEM simulation of an underwater collapse




aspect ratio 'a' of 6 on a slope of 5*

Mesh-based vs Mesh-free techniques

Material Point Method

Porosity in MPM

Material Point Method

MPM slope failure

MPM slope failure: pore pressure changes

Selborne case study (Alonso et al., 2016)

MPM submarine landslide

(Taka., 2012)

Material Point Method: Implementation

  • Generic templatised C++11/14 code
  • 2D/3D MPM Code
  • Elements: Quadrilateral and Hexahedron
  • Shared memory / thread parallelisation
  • BLAS matrix algebra
  • Visualisation VTK
  • Material models: ILE, MC, CamClay, Bingham
  • FE Base library
    • Mesh class
    • Solver class
    • Assembler class
    • MPI parallelisation
    • PETSc matrix algebra / solver
  • FE THM coupling
  • HDF5 support
  • GMSH meshing
  • Two-phase and Two-point formulation

C++ Factory implementation

              
	        void MyElementFactory::RegisterFactoryFunction(
 	        std::string name, std::function classFactoryFunction) {
                  // register the class factory function
                    factoryFunctionRegistry[name] = classFactoryFunction;
                  }
                  std::shared_ptr MyElementFactory::Create(std::string name) {
                    Element *instance = nullptr;

                    // find name in the registry and call factory method.
                    auto it = factoryFunctionRegistry.find(name);
                    if (it != factoryFunctionRegistry.end())
                    instance = it->second();

                    // wrap instance in a shared ptr and return
                    if (instance != nullptr)
                      return std::shared_ptr(instance);
                    else
                      return nullptr;
                  }
              
              

MPM main (snippet)

              
// locate particles in mesh
mesh.locateParticlesAndElements(particleCloudSoilPtr);

// calculate and store shape functions and gradient shape functions
particleCloudSoil.iterateOverParticles(
  boost::bind(&SoilParticleType::evaluateShpfunAndGradShpfunHorMesh, _1));

// calculate nodal acceleration
particleCloudSoil.iterateOverParticles(
  boost::bind(&SoilParticleType::assignForceToNodes, _1));

mesh.iterateOverNodesOfSoilP(
  boost::bind(&NodeType::computeSoilAcceleration, _1, dt, boundaryMiu));
              
              

Iterate over particle

              
//------------------------------------------------------------------------------
// Apply a given functor to all particles making use of std::for_each
// \tparam Toper          Type of the particle operation
// \tparam MPM_ITEMS      MPM Handles
// \tparam PARTICLE_TYPE  Type of particle
// \param[in|out] oper    Specific operation applied to all particles
// \retval Toper          for_each returns the operator

template < typename MPM_ITEMS, typename PARTICLE_TYPE >
template < typename Toper >
Toper mpm::ParticleCloudMpm< MPM_ITEMS, PARTICLE_TYPE >::iterateOverParticles(
    Toper oper) {
  tbb::parallel_for_each(particles_.begin(), particles_.end(),
                         oper);
}
              
              

Particle

              
//------------------------------------------------------------------------------
// \brief class for soil particles
// \details Includes basic calculations for particles
// \tparam MPM_ITEMS   the main template file which includes all the classes

template < typename MPM_ITEMS >
class mpm::ParticleSoilMpm : public mpm::ParticleMpm< MPM_ITEMS > {

  protected:
  using ParticleType::id_;        // identifiying particle number
  using ParticleType::coord_;     // coordinates
  using ParticleType::sfunc_;     // shape functions
  using ParticleType::dsfuncDX_;  // gradient shape function
  using ParticleType::ePtrOfP_;   // pointer to the element which has particle
  using ParticleType::nodesOfParticle_;  // nodes of element in which particle is in

  double mass_;                          //!< mass of particle (constant)
  double density_;           // phase averaged density of soil (=(1-n)density)
  double porosity_;          // porosity of particle
  double phi_;               // friction angle of soil (initial stress calc)
  MaterialPtr material_;     // Pointer to material object
  VecDof velocity_;          // velocity
  VecDof acceleration_;      // acceleration
  Vec6x1 strain_;            // strain
  Vec6x1 dStrain_;           // change of strain in a time step
  Vec6x1 stress_;            // stress
}
              
              

Thank you!!


Krishna Kumar, kks32@cam.ac.uk