Standard.hpp

Go to the documentation of this file.

// This file is part of GranOO, a workbench for DEM simulation.
//
// Author(s)    : - Damien Andre       IRCER/UNILIM, Limoges France
//                  <damien.andre@unilim.fr>
//                - Jean-luc Charles   Arts et Metiers ParisTech, CNRS, I2M, Bordeaux France
//                  <jean-luc.charles@ensam.eu>
//                - Jeremie Girardot   Arts et Metiers ParisTech, CNRS, I2M, Bordeaux France
//                  <jeremie.girardot@ensam.eu>
//                - Cedric Hubert      LAMIH/UPHF, Valenciennes France
//                  <cedric.hubert@uphf.fr>
//                - Ivan Iordanoff     Arts et Metiers ParisTech, CNRS, I2M, Bordeaux France
//                  <ivan.iordanoff@ensam.eu>
//
// Copyright (C) 2008-2019 D. Andre, JL. Charles, J. Girardot, C. Hubert, I. Iordanoff
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 
#ifndef _libDEM_ContactLaw_Standard_hpp_
#define _libDEM_ContactLaw_Standard_hpp_
 
 
#include "GranOO3/Core/ObjectFactory.hpp"
 
#include "GranOO3/DEM/DiscreteElement.hpp"
#include "GranOO3/DEM/SupportShape.hpp"
#include "GranOO3/DEM/Bond.hpp"
#include "GranOO3/DEM/Bond/Spring.hpp"
#include "GranOO3/DEM/ContactLaw.hpp"
#include "GranOO3/DEM/Util/Formula.hpp"
#include "GranOO3/Physic/EnergyBalance.hpp"
#include "GranOO3/FEM/Surface.hpp"
 
#include "GranOO3/Core/SetOf.hpp"
 
#include "GranOO3/Math/Expression.hpp"
 
namespace GranOO3
{
  namespace DEM
  {
 
    template <class T>
    class ContactLaw_Standard : public ContactLaw<T>
    {
 
      GRANOO_OBJECT_FACTORY(ContactLaw_Standard);
 
    public:
      static std::string class_ID() {return "Standard";}
 
    public:
      //CONSTRUCTORS & DESTRUCTORS
      ContactLaw_Standard();
      virtual ~ContactLaw_Standard();
 
      //USEFULL
      virtual void parse_xml();
      virtual void compute_reaction(DEM::DiscreteElement& de1, T& de2, const Geom::Vector& normal, const double& penetration);
      virtual double compute_critical_time_step() const;
      virtual void info(std::ostream &) const;
 
      void set_stiffness(double);
      void set_restitution_coeff(double);
      void set_static_dry_friction_coeff(double);
      void set_dynamic_dry_friction_coeff(double);
      void set_dry_friction_slope(double);
      void set_dry_friction_coeff(double);
      void set_adhesion_force(double);
 
      double get_stiffness() const;
      double get_restitution_coeff() const;
      double get_static_dry_friction_coeff() const;
      double get_dynamic_dry_friction_coeff() const;
      double get_dry_friction_slope() const;
      double get_dry_friction_coeff() const;
      double get_adhesion_force() const;
 
    private:
      ContactLaw_Standard(const ContactLaw_Standard&) = delete;            
      ContactLaw_Standard& operator=(const ContactLaw_Standard&) = delete; 
 
    private:
      bool _exclude_bonded_element;
      double _stiffness;
      double _restitution_coeff;
      double _damping_factor;
      double _static_dry_friction_coeff;
      double _dynamic_dry_friction_coeff;
      double _dry_friction_slope;
      double _dry_friction_coeff;    //old version
      std::string _regularization_type;
      double _dry_friction_life_time;
      double _adhesion_force;
      Math::Expression _normal_force;
      bool _normal_force_mode;
    };
 
    template <> void
    ContactLaw_Standard<FEM::Surface>::compute_reaction(DEM::DiscreteElement &de, FEM::Surface &surface, const Geom::Vector& normal, const double& penetration);
 
    template <class T> void
    ContactLaw_Standard<T>::set_stiffness(double val) {
      UserAssert(val >= 0., "The stiffness must be positive");
      _stiffness = val;
    }
 
    template <class T> void
    ContactLaw_Standard<T>::set_restitution_coeff(double val) {
      UserAssert(val >= 0. && val <= 1.,
                 "The RestitutionCoeff must be between [0,1]");
      _restitution_coeff = val;
      _damping_factor = DEM::restitution_coeff_to_damping_factor(_restitution_coeff);
    }
 
    template <class T> void
    ContactLaw_Standard<T>::set_static_dry_friction_coeff(double val) {
      UserAssert(val >= 0. && val <= 1.,
                 "The StaticDryFrictionCoeff must be between [0,1]");
      _static_dry_friction_coeff = val;
    }
 
    template <class T> void
    ContactLaw_Standard<T>::set_dynamic_dry_friction_coeff(double val) {
      UserAssert(val >= 0. && val <= 1.,
                 "The DynamicDryFrictionCoeff must be between [0,1]");
      _dynamic_dry_friction_coeff = val;
    }
 
    template <class T> void
    ContactLaw_Standard<T>::set_dry_friction_slope(double val) {
      _dry_friction_slope = val;
    }
 
    template <class T> void
    ContactLaw_Standard<T>::set_adhesion_force(double val) {
      _adhesion_force = val;
    }
 
    template <class T> double
    ContactLaw_Standard<T>::get_stiffness() const {
      return _stiffness;
    }
 
    template <class T> double
    ContactLaw_Standard<T>::get_restitution_coeff() const {
      return _restitution_coeff;
    }
 
    template <class T> double
    ContactLaw_Standard<T>::get_static_dry_friction_coeff() const {
      return _dynamic_dry_friction_coeff;
    }
 
    template <class T> double
    ContactLaw_Standard<T>::get_dynamic_dry_friction_coeff() const {
      return _dynamic_dry_friction_coeff;
    }
 
    template <class T> double
    ContactLaw_Standard<T>::get_dry_friction_slope() const {
      return _dry_friction_slope;
    }
 
    template <class T> double
    ContactLaw_Standard<T>::get_adhesion_force() const {
      return _adhesion_force;
    }
 
    template <class T>
    ContactLaw_Standard<T>::ContactLaw_Standard()
      : ContactLaw<T>(),
      _exclude_bonded_element(false),
      _stiffness(0.),
      _restitution_coeff(-1.),
      _damping_factor(-1.),
      _static_dry_friction_coeff(0.),
      _dynamic_dry_friction_coeff(0.),
      _dry_friction_slope(0),
      _dry_friction_coeff(0),
      _regularization_type(""),
      _dry_friction_life_time(0.),
      _adhesion_force(0.),
      _normal_force(),
      _normal_force_mode(false) {
      Core::Mapped<Math::Variable>::Map& m = Math::Variable::get_map();
 
      if (m.count("d") == false) {
        double temp;
        new Math::Variable("d", temp, "InterPenetration");
      }
    }
 
    template <class T>
    ContactLaw_Standard<T>::~ContactLaw_Standard() {
    }
 
    template <class T> void
    ContactLaw_Standard<T>::parse_xml() {
      Core::XmlParser& parser = Core::XmlParser::get();
      parser.read_attribute(Attr::GRANOO_OPTIONAL, "Stiffness", _stiffness);
 
      std::string FnExpr = "";
      parser.read_attribute(Attr::GRANOO_OPTIONAL, "NormalForce", FnExpr);
 
      parser.read_attribute(Attr::GRANOO_OPTIONAL, "RestitutionCoeff", _restitution_coeff);
      parser.read_attribute(Attr::GRANOO_OPTIONAL, "damping_factor", _damping_factor);
      parser.read_attribute(Attr::GRANOO_OPTIONAL, "ExcludeBondedDiscreteElements", _exclude_bonded_element);
      parser.read_attribute(Attr::GRANOO_OPTIONAL, "StaticDryFrictionCoeff", _static_dry_friction_coeff);
      parser.read_attribute(Attr::GRANOO_OPTIONAL, "DynamicDryFrictionCoeff", _dynamic_dry_friction_coeff);
      parser.read_attribute(Attr::GRANOO_OPTIONAL, "DryFrictionSlope", _dry_friction_slope);
      parser.read_attribute(Attr::GRANOO_OPTIONAL, "DryFrictionCoeff", _dry_friction_coeff);
      parser.read_attribute(Attr::GRANOO_OPTIONAL, "RegularizationType", _regularization_type);
      parser.read_attribute(Attr::GRANOO_OPTIONAL, "DryFrictionLifeTime", _dry_friction_life_time);
      parser.read_attribute(Attr::GRANOO_OPTIONAL, "AdhesionForce", _adhesion_force);
 
      // manage the stifness or the normal force function
      XmlAssert((_stiffness==0. && FnExpr!="") || (_stiffness>0. && FnExpr==""),
                "You cannot use both Stiffness and NormalForce !");
 
      if (FnExpr!="") {
        _normal_force_mode = true;
        XmlAssert(_damping_factor==-1.,"Using the NormalForce, you cannot define a damping or restitution factor !");
        XmlAssert(_restitution_coeff==-1.,"Using the NormalForce, you cannot define a damping or restitution factor !");
 
        if (FnExpr.find("d") == std::string::npos)   // not necessary usefull
          FnExpr += "*d";
 
        _normal_force.set(FnExpr);
      }
 
 
      // XmlAssert((_damping_factor == -1 && _restitution_coeff >= 0.) || (_damping_factor >= 0. && _restitution_coeff == -1),
      //           "You cannot use both RestitutionCoeff and damping_factor !");
      XmlAssert(_restitution_coeff == -1 || (_restitution_coeff >= 0. && _restitution_coeff <= 1.),
                "The RestitutionCoeff must be in [0, 1]");
      XmlAssert(_damping_factor    == -1 || (_damping_factor    >= 0. && _damping_factor    <= 1.),
                "The damping_factor must be in [0, 1]");
 
      if (_restitution_coeff >= 0.) {
        // if _restitution_coeff was given, compute _damping_factor:
        _damping_factor = DEM::restitution_coeff_to_damping_factor(_restitution_coeff);
      } else if (_damping_factor >= 0.) {
        // _damping_factor_ was given, compute _restitution_coeff:
        _restitution_coeff = DEM::damping_factor_To_restitution_coeff(_damping_factor);
      }
 
 
      XmlAssert(_regularization_type == "" || _regularization_type == "piecewise" || _regularization_type == "exponential",
                "Available RegularizationTypes are 'piecewise' and 'exponential'");
 
      if (_regularization_type == "") {
        XmlAssert(_dry_friction_coeff >= 0. && _dry_friction_slope >= 0, "A DryFrictionCoeff and a DryFrictionSlope are required.");
        XmlAssert(_dry_friction_coeff >= 0. && _dry_friction_coeff <= 1,  "The DryFrictionCoeff must be in range [0,1].");
        XmlAssert(_dry_friction_slope >= 0, "The DryFrictionSlope must be positive.");
      } else {
        XmlAssert(_static_dry_friction_coeff  >= 0. && _static_dry_friction_coeff  <= 1., "The StaticDryFrictionCoeff must be in range [0,1].");
        XmlAssert(_dynamic_dry_friction_coeff >= 0. && _dynamic_dry_friction_coeff <= 1., "The DynamicDryFrictionCoeff must be in range [0,1].");
        if (_regularization_type == "piecewise")
          XmlAssert(_dry_friction_slope >0, "Piecewise regularization requires a positive DryFrictionSlope");
        if (_regularization_type == "exponential" )
          XmlAssert(_dry_friction_life_time >0, "Exponential regularization requires a positive DryFrictionLifeTime");
 
      }
    }
 
    template <class T> void
    ContactLaw_Standard<T>::info(std::ostream & out) const {
      out << "ContactLaw_Standard"          << granoo::endl;
      out << "\t stiffness              : " << _stiffness << granoo::endl;
      out << "\t restitutionCoeff       : " << _restitution_coeff  << granoo::endl;
      out << "\t dampingFactor          : " << _damping_factor  << granoo::endl;
      out << "\t regularizationType     : " << _regularization_type  << granoo::endl;
      out << "\t staticDryFrictionCoeff : " << _static_dry_friction_coeff  << granoo::endl;
      out << "\t dynamicDryFrictionCoeff: " << _dynamic_dry_friction_coeff  << granoo::endl;
      out << "\t dryFrictionSlope       : " << _dry_friction_slope  << granoo::endl;
      out << "\t dryFrictionCoeff       : " << _dry_friction_coeff  << granoo::endl;
      out << "\t dryFrictionLifeTime    : " << _dry_friction_life_time  << granoo::endl;
      out << "\t adhesionForce          : " << _adhesion_force  << granoo::endl;
    }
 
    template <class T> double
    ContactLaw_Standard<T>::compute_critical_time_step() const {
      // const double minMass = set.lightest().get_mass();
      double minMass = std::numeric_limits<double>::max();
 
      Core::SetOf< DEM::DiscreteElement >& set = Core::SetOf< DEM::DiscreteElement >::get();
      for (unsigned int i = 0; i<set.size(); i++) {
        minMass = set(i).get_mass() < minMass ? set(i).get_mass()  : minMass;
      }
      return sqrt(minMass/get_stiffness());
    }
 
    // contact between 1 DE and a phantom body
    // -----------------------------------------------------------------
 
 
    template <class T> void
    ContactLaw_Standard<T>::compute_reaction(DEM::DiscreteElement& de1,
                                             T& de2,
                                             const Geom::Vector& normal,
                                             const double& penetration) {
 
      if (_exclude_bonded_element) {
        DEM::Element* e1 =  dynamic_cast<DEM::Element*>(&de1);
        DEM::Element* e2 =  dynamic_cast<DEM::Element*>(&de2);
        if (e1->is_bonded(*e2))
          return;
      }
 
      SafeModeAssert(penetration >= 0, "The penetration between de1 and de2 must be positve");
 
      // Compute stiffness reaction
      Geom::Vector Fn;
      if (_normal_force_mode) {
        Math::Variable::get("d").set_val(penetration);
        Fn = normal*_normal_force.value();
      } else
        Fn = normal*_stiffness*penetration;
 
      APPLY_FORCE(Fn, -Fn, "[ContactLaw_Standard]_Elas");
 
      // Compute static and dynamic friction reaction
      if (_static_dry_friction_coeff > 0. || _dynamic_dry_friction_coeff > 0. || _dry_friction_coeff > 0) {
 
        const Geom::Vector AO1   = -de1.get_radius()*normal;
        const Geom::Vector O1O2  =  de1.get_position() - de2.get_position();
        const Geom::Vector AO2   =  AO1+O1O2;
        const Geom::Vector & V1  =  de1.get_linear_velocity();
        const Geom::Vector & V2  =  de2.get_linear_velocity();
        const Geom::Vector & W1  =  de1.get_angular_velocity();
        const Geom::Vector & W2  =  de2.get_angular_velocity();
        Geom::Vector Vrel  = V1 + (AO1^W1) - V2 - (AO2^W2);
        Vrel -= normal*(Vrel*normal);                   // relative velocity in tangential plan
        const Geom::Vector tg = Vrel.unit();       // tangent direction vector
 
        double mu = 0;
        const double VrelModulus = Vrel.norm();
 
        if (_regularization_type == "piecewise") {
          // piecewise linear law for Coulomb Law regularization
          if (VrelModulus < _static_dry_friction_coeff/_dry_friction_slope) {
            mu = _dry_friction_slope*VrelModulus;
          } else if ( (VrelModulus > (_static_dry_friction_coeff/_dry_friction_slope)) && (VrelModulus < ((2*_static_dry_friction_coeff-_dynamic_dry_friction_coeff)/_dry_friction_slope)) ) {
            mu = _static_dry_friction_coeff - _dry_friction_slope*(VrelModulus - _static_dry_friction_coeff/_dry_friction_slope);
          } else if (VrelModulus > ((2*_static_dry_friction_coeff-_dynamic_dry_friction_coeff)/_dry_friction_slope)) {
            mu = _dynamic_dry_friction_coeff;
          }
        } else if (_regularization_type == "exponential") {
          // Decreasing exponential law for Coulomb Law regularization
          mu = _dynamic_dry_friction_coeff + (_static_dry_friction_coeff-_dynamic_dry_friction_coeff)*exp(-VrelModulus/_dry_friction_life_time);
        } else {
          // Old version
          mu = (_dry_friction_slope*VrelModulus < _dry_friction_coeff ) ? _dry_friction_slope*VrelModulus : _dry_friction_coeff;
        }
 
        Geom::Vector Ft = mu * Fn.norm() * tg;
        APPLY_FORCE(-Ft, Ft, "[ContactLaw_Standard]_Fric");
        APPLY_TORQUE(-(-AO1)^Ft, (-AO2)^Ft, "[ContactLaw_Standard]_Fric");
      }
 
      // computes adhesion
      if (_adhesion_force > 0.) {
        const Geom::Vector Fa = _adhesion_force*normal;
        APPLY_FORCE(Fa, -Fa, "[Contactlaw_Standard]_Adhe");
      }
 
      // computes damping reaction
      if (_restitution_coeff != 0.) {
        const double M1 = de1.get_mass();
        const double M2 = de2.get_mass();
        const double M  = M1*M2/(M1+M2);
        const double K  = _stiffness;
        const double dampCoeff     = 2.*_damping_factor*sqrt(K*M);
        const double velocityDiff  = (de2.get_linear_velocity() - de1.get_linear_velocity()) * normal;
        const Geom::Vector Fd = dampCoeff * velocityDiff * normal;
 
        APPLY_FORCE(Fd, -Fd, "[ContactLaw_Standard]_Damp");
      }
    }
 
    template<> inline void
    ContactLaw_Standard<DEM::SupportShape>::compute_reaction(DEM::DiscreteElement& de,
                                                             DEM::SupportShape& de2,
                                                             const Geom::Vector& normal,
                                                             const double& penetration) {
      SafeModeAssert(penetration >= 0, "The penetration between de1 and de2 must be positve");
      Physic::EnergyBalance& energyBalance = Physic::EnergyBalance::get();
 
      // Compute stiffness reaction
      Geom::Vector Fn;
      if (_normal_force_mode) {
        Math::Variable::get("d").set_val(penetration);
        Fn = normal*_normal_force.value();
      } else
        Fn = normal*_stiffness*penetration;
 
      de.force() += Fn;
      if (energyBalance.is_enable())
        energyBalance.add_force_on(de, "[ContactLaw_Standard]_Elas", Fn);
 
      // Compute static and dynamic friction reaction
      if (_static_dry_friction_coeff > 0. || _dynamic_dry_friction_coeff > 0. || _dry_friction_coeff > 0) {
        const Geom::Vector AO1   = -de.get_radius()*normal;
        const Geom::Vector O1O2  =  de.get_position() - de2.get_position();
        const Geom::Vector AO2   =  AO1+O1O2;
        const Geom::Vector & V1  =  de.get_linear_velocity();
        const Geom::Vector & V2  =  de2.get_linear_velocity();
        const Geom::Vector & W1  =  de.get_angular_velocity();
        const Geom::Vector & W2  =  de2.get_angular_velocity();
 
        Geom::Vector Vrel  = V1 + (AO1^W1) - V2 - (AO2^W2);
        Vrel -= normal*(Vrel*normal);                   // relative velocity in tangential plan
 
        double mu = 0;
        const double VrelModulus = Vrel.norm();
 
        if (_regularization_type == "piecewise") {
          // piecewise linear law for Coulomb Law regularization
          if (VrelModulus < _static_dry_friction_coeff/_dry_friction_slope) {
            mu = _dry_friction_slope*VrelModulus;
          } else if ( (VrelModulus > (_static_dry_friction_coeff/_dry_friction_slope)) &&
                      (VrelModulus < ((2*_static_dry_friction_coeff-_dynamic_dry_friction_coeff)/_dry_friction_slope)) ) {
            mu = _static_dry_friction_coeff - _dry_friction_slope*(VrelModulus - _static_dry_friction_coeff/_dry_friction_slope);
          } else if (VrelModulus > ((2*_static_dry_friction_coeff-_dynamic_dry_friction_coeff)/_dry_friction_slope)) {
            mu = _dynamic_dry_friction_coeff;
          }
        } else if (_regularization_type == "exponential") {
          // Decreasing exponential law for Coulomb Law regularization
          mu = _dynamic_dry_friction_coeff + (_static_dry_friction_coeff-_dynamic_dry_friction_coeff)*exp(-VrelModulus*1e8);
        } else {
          // Old version
          mu = (_dry_friction_slope*VrelModulus < _dry_friction_coeff ) ? _dry_friction_slope*VrelModulus : _dry_friction_coeff;
        }
 
        Geom::Vector Ft =  mu * Fn.norm() * Vrel.unit();
 
        de.force()  -= Ft;
        de.torque() -= (-AO1)^Ft;
 
        if (energyBalance.is_enable()) {
          energyBalance.add_force_on (de, "[ContactLaw_Standard]_Fric", -Ft);
          energyBalance.add_torque_on(de, "[ContactLaw_Standard]_Fric", AO1^Ft);
        }
      }
 
      // computes adhesion
      if (_adhesion_force > 0.) {
        const Geom::Vector Fa = _adhesion_force*normal;
        de.force() += Fa;
 
        if (energyBalance.is_enable())
          energyBalance.add_force_on(de, "[ContactLaw_Standard]_Adhe", Fa);
      }
 
      // computes damping reaction
      if (_restitution_coeff != 0.) {
        const double M = de.get_mass();
        const double K  = _stiffness;
        const double dampCoeff     = 2.*_damping_factor*sqrt(K*M);
        const double velocityDiff  =  -de.get_linear_velocity() * normal;
        const Geom::Vector Fd = dampCoeff * velocityDiff * normal;
        de.force() += Fd;
 
        if (energyBalance.is_enable())
          energyBalance.add_force_on(de, "[ContactLaw_Standard]_Damp", Fd);
      }
    }
 
  }
}
 
#endif