Contact Algorithm

[thiagordonho]

Contact Algorithm

Summary

Implementation of a contact algorithm to handle the frictional interaction between adjacent materials at the region of contact between the two, also known as the multi-material contact or interface. This implementation is based on the procedures provided by Bardenhagen et al. (2000) and Nairn (2013).

Motivation

Currently, the CB-Geo MPM code computes a single velocity field for one-phase problems. This means that, upon computing the nodal velocities towards the end of each time step, the code will update the material point’s (or particle’s) position using the same velocity field, regardless of the distinct bodies being represented. The distinct bodies will thus behave as if they are "stuck" to each other. Nevertheless, this is not always the case for every model, since one may desire to attribute to the body the capacity of moving freely if separating, or of developing friction at the contact surface with another body if they are approaching.

The proposed feature, thus, allows one to develop simulations where multi-material interface interaction is essential in the analysis. Cases such as soil-structure interaction, slip surfaces in slope stability analyses, fracture movement/rupture, among others, represent common geotechnical problems that fit this contact requirement. By applying this feature, the user is able to assign the desired friction characteristics between every pair of bodies.

Design Detail

The MPM contact algorithm as described by Bardenhagen et al. (2000) can be summarized into the following stages: reset the nodal properties, map the particle’s properties to the nodes, compute the nodal kinematics of separate bodies by satisfying the conservation of linear momentum, correct the kinematics by applying contact law and no-interpenetration conditions at the interface between bodies, and update the particles’ kinematics. This brief overview of the contact algorithm stages also indicates that the nodal properties must be stored for separate bodies, instead of the singular nodal velocity field of the one-phase conventional MPM. Therefore, a data structure solution is included in this design detail to explain the chosen approach for implementing the algorithm to the current CB-Geo MPM code. The stages of the representative contact time-step loop are fulfilled by iterating over the particles or nodes to perform each of the tasks. Hence, new methods, hereafter referred to as contact-related methods, are added to the Particle, Node, and Mesh classes to perform these tasks. An additional Contact class is proposed to store the parameters of the contact law and to orderly call the contact-related methods. Finally, a description of the input needed for using this feature is given.

Data structure

As previously stated, the stored nodal properties in this scheme must differentiate between distinct bodies. Therefore, the suggested approach is to create a NodalProperty struct representing a property pool that will be responsible for storing and manipulating the data of each nodal property in the contact approach. The pool is declared and handled within the Mesh class under the alias nodal_properties

//! Nodal property pool
std::shared_ptr<mpm::NodalProperties> nodal_properties_{nullptr};

In this entity, different properties are handled with a map container, where each property (key) is associated with its own Eigen::Matrix (value). The total number of columns in these matrices is the same as the number of materials in the model. The total number of rows is equal to the number of nodes multiplied by the dimension of the property – 1 if scalar and Tdim if vector.

std::map<std::string, Eigen::MatrixXd> properties_;

Accessing and modifying the values in these matrices is done with implemented methods in the NodalProperties struct. These methods include the property name, node id, and material id as arguments for returning, assigning, and updating property values in the specific property pool. For instance, in a 2D environment consisting of 2 nodes and 3 materials, each vectorial property will be stored in the following structure:

Eigen::Matrix<double,4,3> data2;
data2 <<  p_000, p_001, p_002,
         p_100, p_101, p_102,
         p_010, p_011, p_012,
         p_110, p_111, p_112;

where p_ijk is the momentum in the i direction, at node j and material k. Accessing the property from node 0 and material 1 means that p_001 and p_101 will be accessed.

Before starting the time step iteration of the MPM solver, the property handle is initialized and the key-value pairs are created. Additionally, a handle to nodal_properties_ is also created in each node to avoid passing the created struct object as an argument when its access is necessary from Node objects. This procedure is done within the create_nodal_property() method in Mesh described below.

// Create the nodal properties' map
template <unsigned Tdim>
void mpm::Mesh<Tdim>::create_nodal_properties() {
  // Initialise the shared pointer to nodal properties
  nodal_properties_ = std::make_shared<mpm::NodalProperties>();

  // Check if nodes_ and materials_is empty and throw runtime error if they are
  if (nodes_.size() != 0 && materials_.size() != 0) {
    // Compute number of rows in nodal properties for vector entities
    const unsigned nrows = nodes_.size() * Tdim;
    // Create pool data for each property in the nodal properties struct
    // object. Properties are named in the plural form
    nodal_properties_->create_property("masses", nodes_.size(),
                                       materials_.size());
    nodal_properties_->create_property("momenta", nrows, materials_.size());
    // […] Proceed similarly to create other pairs of the property map

    // Iterate over all nodes to initialise the property handle in each node
    // and assign its node id as the prop id in the nodal property data pool
    for (auto nitr = nodes_.cbegin(); nitr != nodes_.cend(); ++nitr)
      (*nitr)->initialise_property_handle((*nitr)->id(), nodal_properties_);
  } else {
    throw std::runtime_error("Number of nodes or number of materials is zero");
  }
}

The following are the keys (properties) stored in the map object in nodal_properties:

• "masses"
• "momenta"
• "displacements"
• ”domains”
• ”domain_gradients”
• ”normal_unit_vectors”
• ”internal_forces”
• ”external_forces”
• ”current_velocities”
• "velocities"
• ”accelerations”
• ”relative_velocities”

All keys are saved as std::string, representing the name of the property in plural form.

Contact-related methods

The multi-material contact approach proposed by Bardenhagen et al. (2000) can be subdivided into a series of steps in the code. Except for the creation of nodal_properties and each of its key-value pairs, which happens once at the start of the simulation, all steps happen sequentially within the time-step loop:

1. Initialise nodal_properties, setting all the parameters to zero.
2. Initialise material ids set in Node Class by iterating over all nodes. This set of material ids is used in future iterations within the node to compute properties only for the mapped materials, thus accounting for the sparsity in the Eigen Matrices of nodal_properties.
3. Append material ids to the designated set in the Node Class by iterating over all nodes.
4. Map multi-material mass and momentum to the nodal_properties by iterating over all material points.
5. Compute the velocity by iterating over all nodes.
6. Map domain (volume) and gradient of the domain to the nodal_properties by iterating over all material points.
7. Compute the normal unit vector by iterating over all nodes.
8. Map the body forces and internal forces by iterating over all material points.
9. Apply traction on particles and map it to the nodes.
10. Apply concentrated force to all nodes.
11. Compute the acceleration and velocity by satisfying the conservation of mass equation at all nodes but for separate bodies.
12. Compute the relative velocity at all the nodes.
13. Verify if bodies are approaching at contact and apply the contact mechanics if true.
14. Apply particle velocity constraints.
15. Update particle’s kinematics from the derived nodal acceleration and velocities for each respective material.

The steps itemized above will either take place in the particle or node level. They are thus implemented as methods within the Particle or Node classes and are called by iterating over all particles and nodes, respectively. Most of the functions are similar to the conventional MPM but an additional iteration over the material ids is done to account for distinct bodies and the obtained values are stored in nodal_properties. For instance, the function to map mass and momentum is as follows:

//! Map multimaterial properties to nodes
template <unsigned Tdim>
void mpm::Particle<Tdim>::map_multimaterial_mass_momentum_to_nodes() noexcept {
  // Check if particle mass is set
  assert(mass_ != std::numeric_limits<double>::max());

  // Unit 1x1 Eigen matrix to be used with scalar quantities
  Eigen::Matrix<double, 1, 1> nodal_mass;

  // Map mass and momentum to nodal property taking into account the material id
  for (unsigned i = 0; i < nodes_.size(); ++i) {
    nodal_mass(0, 0) = mass_ * shapefn_[i];
    nodes_[i]->update_property(true, "masses", nodal_mass, this->material_id(),
                               1);
    nodes_[i]->update_property(true, "momenta", velocity_ * nodal_mass,
                               this->material_id(), Tdim);
  }
}

The steps where the methods are similar to the original MPM code are steps 4, 5, 8-11. The remaining methods are either entirely new to the code or differ significantly from the corresponding conventional MPM methods.

Steps 1-3 correspond to the initialization of the nodal_properties and material ids that are mapped to the node. They represent the first steps taken within the time-step loop before any computation takes place. Step 1 sets all nodal properties to zero, step 2 clears the set of material ids in every node (happens within initialise() at the Node class) and step 3 has a function of its own to append new material ids to the node as the code iterate through all particles:

//! Assign material id of this particle to nodes
template <unsigned Tdim>
void mpm::Particle<Tdim>::append_material_id_to_nodes() const {
  for (unsigned i = 0; i < nodes_.size(); ++i)
    nodes_[i]->append_material_id(this->material_id());
}

Steps 6 and 7 are the first to introduce a novelty to the code that is essential to the contact algorithm. These steps pertain to the computation of the normal unit vector at every node and with respect to each body. The normal unit vector computation follows the possible cases presented by Nairn (2013). The cases used in the context of this code are the pure volume gradient (PVG), the average volume gradient (AVG), or the maximum volume gradient (MVG). The three methods use the mapped volume (domain) and the volume gradient (domain gradient) that is determined in step 6 similarly to the way other properties are mapped to the nodes. Afterward, the normal is computed in step 7 in one of the three ways:

• PVG: normalizing the mapped domain gradient.
• MVG: if the only two bodies are in contact at the node, normalize the greatest mapped domain gradient between the two and assign its opposite to the other body.
• AVG: the weighted average of the domain gradients with respect to the mapped gradient is taken at the nodes. The normal unit vector is the normalized resulting vector.

The normal computation method used must be user-defined as indicated in the input section of this document.

The core difference of the implementation of the contact algorithm to the code is the application of contact mechanics in steps 12 and 13. Step 12 simply computes the relative velocity at the node and body, which is the difference between the body's nodal velocity and the nodal velocity of the center of mass – the latter is the already computed nodal velocity as in the conventional MPM.

//! Compute multimaterial relative velocities
template <unsigned Tdim, unsigned Tdof, unsigned Tnphases>
void mpm::Node<Tdim, Tdof,
               Tnphases>::compute_multimaterial_relative_velocity() {
  // Iterate over all materials in the material_ids set
  node_mutex_.lock();
  for (auto mitr = material_ids_.begin(); mitr != material_ids_.end(); ++mitr) {
    VectorDim velocity_material =
        property_handle_->property("velocities", prop_id_, *mitr, Tdim);
    VectorDim relative_velocity = velocity_material - this->velocity_.col(0);

    // Assign relative velocities
    property_handle_->assign_property("relative_velocities", prop_id_, *mitr,
                                      relative_velocity, Tdim);
  }
  node_mutex_.unlock();
}

In Step 13, the contact-mechanics as described in Bardenhagen et al. (2000) is applied. At first, the approaching condition is verified. If the relative velocity’s normal component is greater than 0, the body is approaching another. Otherwise, it is considered to be separating. Only the former condition will result in corrections to the body’s nodal velocity.

      double velocity_normal = relative_velocity.dot(normal_unit_vector);

      // Check if the material is approaching the other materials (v_norm > 0)
      if (velocity_normal > 0) {

The first correction represents the no-interpenetration of the bodies, and it translates as the deduction of the normal component of the relative velocity from the body’s current nodal velocity.

        // Compute the normal correction
        VectorDim normal_correction = -velocity_normal * normal_unit_vector;

A second correction is applied to account for the frictional relationship between the bodies. This correction checks if the tangential component of the relative velocity exceeds the maximum frictional deduction at the node. If it does, this maximum friction is deducted from the body’s current nodal velocity. Otherwise, the tangential component of the relative velocity is deducted instead. This procedure was simply equated as shown in the figure below. An if-else statement was included in the code to deal with the tangential corrections in 2D or 3D.

Finally, the corrected nodal acceleration is assigned as the difference between the corrected nodal velocity and the initial nodal velocity, all divided by the time step. Moreover, the initial nodal velocity corresponds to the nodal velocity derived from mapping the particles' mass and momentum to the nodes (steps 4 and 5).

Contact class

To store the specific contact-related parameters, a Contact class is used in this implementation. This class also contains the ordered calls to the contact-related methods. This facilitates the future implementation of other contact laws (not just frictional), minimizing the changes to ideas and structure of the contact algorithm. The proceedings mentioned so far are all part of the ContactFriction class, derived from Contact. A pointer to the contact object is declared in the MPMBase class. It is initialized and called within the derived MPMExplicit::solver() method according to the input JSON file. The contact-related methods are used in separate methods of the ContactFriction to permit future changes to specific sections and to fit the order of calls to these functions with the already existing ones in the conventional MPM scheme. The following is a code snippet of the ContactFriction class:

//! ContactFriction class
//! \brief ContactFriction base class
//! \tparam Tdim Dimension
template <unsigned Tdim>
class ContactFriction : public Contact<Tdim> {
 public:
  //! Default constructor with mesh class
  ContactFriction(const std::shared_ptr<mpm::Mesh<Tdim>>& mesh, double friction,
                  std::string normal_type);

  //! Intialize
  virtual inline void initialise() override;

  //! Compute nodal kinematics
  virtual inline void compute_nodal_kinematics() override;

  //! Compute contact forces
  //! \param[in] gravity Gravity vector
  //! \param[in] phase Index to indicate material phase
  //! \param[in] time Current time in the simulation
  //! \param[in] concentrated_nodal_forces Boolean for if a concentrated force
  //! is applied or not
  virtual inline void compute_contact_forces(
      const Eigen::Matrix<double, Tdim, 1>& gravity, unsigned phase,
      double time, bool concentrated_nodal_forces) override;

  //! Compute contact nodal kinematics
  //! \param[in] dt Timestep in analysis
  virtual inline void compute_contact_kinematics(double dt) override;

  //! Update particle position
  //! \param[in] dt Timestep in analysis
  //! \param[in] velocity_update Update particle velocity from nodal vel
  virtual inline void update_particles_contact(double dt,
                                               bool velocity_update) override;

 protected:
  //! Mesh object
  using mpm::Contact<Tdim>::mesh_;
  //! Coefficient of friction
  double friction_{.0};
  //! Type of normal detection
  std::string normal_type_{"default"};
};  // Contactfriction class

Input

To consider frictional contact in a model, the user must define the friction coefficient and the type of normal computation within the analysis object of the JSON file. For example, a coefficient of friction of 0.2 and the normal computation type of MVG are to be used, the JSON file would include the following:

"analysis": {
    "interface": {
        "interface": true,
        "friction": 0.2,
        "normal_type": "MVG"
    }
}

Notice that the use of contact can be toggled on and off by setting the "interface" entry to true or false, respectively.

Additionally, the user must define a new material with an id for each body that they intend to use in the model. The material ids used must start at 0 and must also be consecutive. The material parameters can be the same between materials, as long as different ids are used to represent distinct bodies.

Currently, the proposed code does not allow for the use of multiple coefficients of friction between bodies. Regardless of the material ids of the bodies that are approaching, the code will apply the same coefficient of friction and normal type detection.

Drawbacks

• It will be more expensive computationally.
• It will require a lot of memory to store the information for every node, especially for simulations that comprise a computational grid with many nodes.
• Memory stored might not be used in the case of nodes in the mesh that have one or no materials adjacent to it.
• The code currently only allows for the use of one contact law (frictional) and one set of parameters (coefficient of friction and normal computation type).

Rationale and Alternatives

Another alternative to the Algorithm implemented would be to directly store the information in the nodes instead of having a property pool (nodal_properties). This can even encompass dynamic allocation of memory to avoid the unnecessary usage of memory for the nodes with one or no adjacent materials. However, this alternative would result in a slower algorithm with the possibility of memory leaks happening if no proper release of memory is considered.

Prior Art

This implementation is highly based on the directions provided by Bardenhagen et al. (2001) and Nairn (2013). In fact, the former reference is the one that introduces the portion of the contact algorithm envisioned for this design, whereas the latter introduced the Imperfect Interface Theory to deal with the cases where adjacent materials are not in contact but still constitute what the author calls imperfect interface. However, the details about the imperfect interface contribution is out of the scope of this design. The latter also presents additional directions on computing the normal unit vectors.

Unresolved questions

• A proper algorithm for identifying the interfaces may have to be implemented in order to drastically reduce the amount of memory used and update the size of nodal_properties to account only for the nodes where more than one material is identified to be adjacent to the node.

[kks32]

Yes, I'm happy with this, I move to vote. 👍 or 💬 if you like to discuss more.

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

[kks32]

We are currently waiting for benchmarks and validation problems.

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

[kks32]

@thiagordonho

[stale]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.