yade-dev team mailing list archive

[Behzad Majidi <behzad.majidi@xxxxxxxxx>]

Hi there,

I have developed a new contact model, namely Burger's viscoelastic model
for YADE. I have tried and verified the model in my own YADE trunk.
I'd like to publish it on the main trunk to get it reviewed and further
developed (if any)

 was not able to commit the changes directly on github, using git push
command. My username is bmajidi


However, I have attached the model files to this email and I appreciate if
you can commit the changes for me.

Thank you,


-- 
Behzad Majidi

Research Student
Aluminium Research Centre
REGAL-Université Laval
Québec, QC, Canada





Please consider the environment before printing this E-mail.




-- 
Behzad Majidi

Research Student
Aluminium Research Centre
REGAL-Université Laval
Québec, QC, Canada





Please consider the environment before printing this E-mail.

// 2014 © Behzad Majidi <behzad.majidi@xxxxxxxxx>
// 2014 © Ricardo Pieralisi <ricpieralisi@xxxxxxxxx>

#include"Cohburgers.hpp"
#include<core/State.hpp>
#include<core/Interaction.hpp>
#include<pkg/dem/ScGeom.hpp>
#include<core/Omega.hpp>
#include<core/Scene.hpp>
#include<pkg/dem/Shop.hpp>
#include<pkg/common/InteractionLoop.hpp>
#include<pkg/common/Sphere.hpp>
#include<pkg/common/Facet.hpp>
#include<pkg/common/Wall.hpp>


YADE_PLUGIN((CohBurgersMat)(CohBurgersPhys)(Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys)(Ip2_CohFrictMat_CohBurgersMat_CohBurgersPhys)(Ip2_FrictMat_CohBurgersMat_CohBurgersPhys)(Law2_ScGeom6D_CohBurgersPhys_CohesiveBurgers)); 



//********************** Ip2_FrictMat_CohBurgersMat_CohBurgersPhys ****************************/

void Ip2_FrictMat_CohBurgersMat_CohBurgersPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) 
{

  if(interaction->phys) return;
  shared_ptr<CohBurgersPhys> phys (new CohBurgersPhys());
 

  int i1 = b1->getClassIndex(); // get actual material index of b1
  int i2 = b2->getClassIndex(); // get actual material index of b2
  int cbmi = CohBurgersMat::getClassIndexStatic(); // get index of CohBurgersMat
  FrictMat* mat1;
  CohBurgersMat* mat2;
  if (i2 == cbmi) { // b2 is CohburgersMat and b1 is FrictMat
  mat1 = dynamic_cast<FrictMat*>(b1.get());
  mat2 = dynamic_cast<CohBurgersMat*>(b2.get());
  } else if (i1 == cbmi) { // b1 is CohburgersMat and b2 is FrictMat
  mat1 = dynamic_cast<FrictMat*>(b2.get());
  mat2 = dynamic_cast<CohBurgersMat*>(b1.get());
  } else { // should not happen, but to be sure..
  LOG_FATAL("TODO");
  }


  ScGeom6D* geom= YADE_CAST<ScGeom6D*>(interaction->geom.get());

  Real D1 = geom->radius1*2;
  Real D2 = geom->radius2*2;

//Load the rheological parameters
  
  Real kn1 = (mat1->young)*D1; Real ks1 = (mat1->young)*D1*(mat1->poisson);
  Real kmn2 = (mat2->kmn)*D2; Real kkn2 = (mat2->kkn)*D2;
  Real cmn2 = (mat2->cmn)*D2; Real ckn2 = (mat2->ckn)*D2;
  Real kms2 = (mat2->kms)*D2; Real kks2 = (mat2->kks)*D2;
  Real cms2 = (mat2->cms)*D2; Real cks2 = (mat2->cks)*D2;

  phys->kmn = kn1;
  phys->kms = ks1;
  phys->cmn = cmn2;
  phys->cms = cms2;

  phys->kkn = kkn2;
  phys->kks = kks2;
  phys->ckn = ckn2;
  phys->cks = cks2;
  
  
  if ((setCohesionOnNewContacts || setCohesionNow) && mat2->isCohesive)
  {
   phys->MaxDef = 0.125*(D1+D2);
   phys->normalAdhesion= (mat2->normalCohesion)*0.5*pow(std::min(D1*0.5,D2*0.5),2);
   phys->shearAdhesion= (mat2->shearCohesion)*0.5*pow(std::min(D1*0.5,D2*0.5),2);
   phys->fragile=mat2->fragile;
   phys->cohesionBroken = false;
   phys->initCohesion   = true;
  }
  interaction->phys = shared_ptr<CohBurgersPhys>(phys);
}



//********************** Ip2_CohFrictMat_CohBurgersMat_CohBurgersPhys ****************************/

void Ip2_CohFrictMat_CohBurgersMat_CohBurgersPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) 
{

  if(interaction->phys) return;
  shared_ptr<CohBurgersPhys> phys (new CohBurgersPhys());
 

  int i1 = b1->getClassIndex(); // get actual material index of b1
  int i2 = b2->getClassIndex(); // get actual material index of b2
  int cbmi = CohBurgersMat::getClassIndexStatic(); // get index of CohBurgersMat
  CohFrictMat* mat1;
  CohBurgersMat* mat2;
  if (i2 == cbmi) { // b2 is CohburgersMat and b1 is FrictMat
  mat1 = dynamic_cast<CohFrictMat*>(b1.get());
  mat2 = dynamic_cast<CohBurgersMat*>(b2.get());
  } else if (i1 == cbmi) { // b1 is CohburgersMat and b2 is FrictMat
  mat1 = dynamic_cast<CohFrictMat*>(b2.get());
  mat2 = dynamic_cast<CohBurgersMat*>(b1.get());
  } else { // should not happen, but to be sure..
  LOG_FATAL("TODO");
  }


  ScGeom6D* geom= YADE_CAST<ScGeom6D*>(interaction->geom.get());

  Real D1 = geom->radius1*2;
  Real D2 = geom->radius2*2;

//Load the rheological parameters
  
  Real kn1 = (mat1->young)*D1; Real ks1 = (mat1->young)*D1*(mat1->poisson);
  Real kmn2 = (mat2->kmn)*D2; Real kkn2 = (mat2->kkn)*D2;
  Real cmn2 = (mat2->cmn)*D2; Real ckn2 = (mat2->ckn)*D2;
  Real kms2 = (mat2->kms)*D2; Real kks2 = (mat2->kks)*D2;
  Real cms2 = (mat2->cms)*D2; Real cks2 = (mat2->cks)*D2;


  phys->kmn = kn1;
  phys->kms = ks1;
  phys->cmn = cmn2;
  phys->cms = cms2;

  phys->kkn = kkn2;
  phys->kks = kks2;
  phys->ckn = ckn2;
  phys->cks = cks2;
  
  
  if ((setCohesionOnNewContacts || setCohesionNow) && mat1->isCohesive && mat2->isCohesive)
  {
   phys->MaxDef = 0.125*(D1+D2);
   phys->normalAdhesion= (mat2->normalCohesion)*0.5*pow(std::min(D1*0.5,D2*0.5),2);
   phys->shearAdhesion= (mat2->shearCohesion)*0.5*pow(std::min(D1*0.5,D2*0.5),2);
   phys->fragile=mat2->fragile;
   phys->cohesionBroken = false;
   phys->initCohesion   = true;
  }
  interaction->phys = shared_ptr<CohBurgersPhys>(phys);
}


//********************** Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys ****************************/

void Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) 
{
  if(interaction->phys) return;
  shared_ptr<CohBurgersPhys> phys (new CohBurgersPhys());

  CohBurgersMat* mat1 = YADE_CAST<CohBurgersMat*>(b1.get());
  CohBurgersMat* mat2 = YADE_CAST<CohBurgersMat*>(b2.get());

  ScGeom6D* geom= YADE_CAST<ScGeom6D*>(interaction->geom.get());

  Real D1 = geom->radius1*2;
  Real D2 = geom->radius2*2;

//Load the rheological parameters

  Real kmn1 = (mat1->kmn)*2*D1; Real kkn1 = (mat1->kkn)*2*D1;
  Real cmn1 = (mat1->cmn)*2*D1; Real ckn1 = (mat1->ckn)*2*D1;
  Real kms1 = (mat1->kms)*2*D1; Real kks1 = (mat1->kks)*2*D1;
  Real cms1 = (mat1->cms)*2*D1; Real cks1 = (mat1->cks)*2*D1;
  Real kmn2 = (mat2->kmn)*2*D2; Real kkn2 = (mat2->kkn)*2*D2;
  Real cmn2 = (mat2->cmn)*2*D2; Real ckn2 = (mat2->ckn)*2*D2;
  Real kms2 = (mat2->kms)*2*D2; Real kks2 = (mat2->kks)*2*D2;
  Real cms2 = (mat2->cms)*2*D2; Real cks2 = (mat2->cks)*2*D2;

  phys->kmn = contactParameterCalculation(kmn1,kmn2);
  phys->kms = contactParameterCalculation(kms1,kms2);
  phys->cmn = contactParameterCalculation(cmn1,cmn2);
  phys->cms = contactParameterCalculation(cms1,cms2);

  phys->kkn = contactParameterCalculation(kkn1,kkn2);
  phys->kks = contactParameterCalculation(kks1,kks2);
  phys->ckn = contactParameterCalculation(ckn1,ckn2);
  phys->cks = contactParameterCalculation(cks1,cks2);

  if ((setCohesionOnNewContacts || setCohesionNow) && mat1->isCohesive && mat2->isCohesive)
  {
   phys->MaxDef = 0.25*(D1+D2);
   phys->normalAdhesion= std::min(mat1->normalCohesion, mat2->normalCohesion)*pow(std::min(D1*0.5,D2*0.5),2);
   phys->shearAdhesion= (mat2->shearCohesion)*0.5*pow(std::min(D1*0.5,D2*0.5),2);
   phys->fragile=(mat1->fragile || mat2->fragile);
   phys->cohesionBroken = false;
   phys->initCohesion   = true;
  }
  interaction->phys = shared_ptr<CohBurgersPhys>(phys);

}

//********************** CohesiveBurgersContactLaw ****************************/

bool Law2_ScGeom6D_CohBurgersPhys_CohesiveBurgers::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I) 
{
  Vector3r force = Vector3r::Zero();
  Vector3r torque1 = Vector3r::Zero();
  Vector3r torque2 = Vector3r::Zero();

  CohBurgersPhys& phys=*YADE_CAST<CohBurgersPhys*>(_phys.get());

  if (computeForceTorqueCohBurgers(_geom, _phys, I, force, torque1, torque2) and (I->isActive)){
    const int id1 = I->getId1();
    const int id2 = I->getId2();
		
    addForce (id1,-force,scene);
    addForce (id2, force,scene);
    addTorque(id1, torque1,scene);
    addTorque(id2, torque2,scene);
    return true;
   } else return false;
}


 
bool computeForceTorqueCohBurgers(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r & torque2) {



  const ScGeom6D& geom=*YADE_CAST<ScGeom6D*>(_geom.get());
  CohBurgersPhys& phys=*YADE_CAST<CohBurgersPhys*>(_phys.get());

  Scene* scene=Omega::instance().getScene().get();

  const int id1 = I->getId1();
  const int id2 = I->getId2();

  const BodyContainer& bodies = *scene->bodies;
    
  const State& de1 = *YADE_CAST<State*>(bodies[id1]->state.get());
  const State& de2 = *YADE_CAST<State*>(bodies[id2]->state.get());

  const Real& dt = scene->dt;
 
  if (phys.fragile && phys.normalAdhesion > 0 && phys.normalForce.norm() > phys.normalAdhesion) {
		phys.cohesionBroken= true;
		return false;	
  }

  if (phys.fragile && std::abs(geom.penetrationDepth) > phys.MaxDef) {
		phys.cohesionBroken= true;
		return false;	
  }

  if (I->isFresh(scene)) {
      phys.ukn = Vector3r(0,0,0);
      phys.uks = Vector3r(0,0,0);
      phys.shearForce = Vector3r(0,0,0);
      phys.normalForce = phys.kmn * geom.penetrationDepth * geom.normal;
      return true;
    } else {
      Real An = 1.0 + phys.kkn * dt / (2.0 * phys.ckn);
      Real Bn = 1.0 - phys.kkn * dt / (2.0 * phys.ckn);
      Real Cn = dt / (2.0 * phys.ckn * An) + 1.0 / phys.kmn + dt / (2.0 * phys.cmn);
      Real Dn = dt / (2.0 * phys.ckn * An) - 1.0 / phys.kmn + dt / (2.0 * phys.cmn);

      Real As = 1.0 + phys.kks * dt / (2.0 * phys.cks);
      Real Bs = 1.0 -  phys.kks * dt / (2.0 * phys.cks);
      Real Cs = dt / (2.0 * phys.cks * As) + 1.0 / phys.kms + dt / (2.0 * phys.cms);
      Real Ds = dt / (2.0 * phys.cks * As) - 1.0 / phys.kms + dt / (2.0 * phys.cms);

      const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero(); 
      const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(I->cellDist): Vector3r::Zero(); 

      const Vector3r c1x = (geom.contactPoint - de1.pos);
      const Vector3r c2x = (geom.contactPoint - de2.pos - shift2);

      const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) + shiftVel;
      const Vector3r normalVelocity = geom.normal.dot(relativeVelocity) * geom.normal;
      const Vector3r shearVelocity = relativeVelocity-normalVelocity;

      const Vector3r normalForceT = (normalVelocity * dt + phys.ukn * (1.0 - Bn / An) - phys.normalForce * Dn) / Cn;
      phys.ukn = (phys.ukn * Bn + (normalForceT + phys.normalForce) * dt / (2.0 * phys.ckn)) / An;
      phys.normalForce = normalForceT;

      const Vector3r shearForceT = -1 * (shearVelocity * dt + phys.uks * (1.0 - Bs / As) - phys.shearForce * Ds) / Cs;
      phys.uks = (phys.uks * Bs + (shearForceT + phys.shearForce) * dt / (2.0 * phys.cks)) / As;
      phys.shearForce = shearForceT;

      const Real maxFs = phys.normalForce.squaredNorm() * std::pow(phys.tangensOfFrictionAngle,2);
      if( phys.shearForce.squaredNorm() > maxFs ) {
		  
        const Real ratio = sqrt(maxFs) / phys.shearForce.norm();
        phys.shearForce *= ratio;
      }

      force = phys.normalForce + phys.shearForce;
      torque1 = -c1x.cross(force);
      torque2 =  c2x.cross(force);
      return true; 
    
  }

}












// 2014 © Behzad Majidi <behzad.majidi@xxxxxxxxx>
// 2014 © Ricardo Pieralisi <ricpieralisi@xxxxxxxxx>

#pragma once

#include<core/Material.hpp>
#include<pkg/dem/FrictPhys.hpp>
#include<pkg/dem/ViscoelasticPM.hpp>
#include<pkg/dem/CohesiveFrictionalContactLaw.hpp>
#include<pkg/common/Dispatching.hpp>
#include<pkg/dem/ScGeom.hpp>
#include<pkg/dem/DemXDofGeom.hpp>
#include<pkg/common/MatchMaker.hpp>


//********************** CohBurgersMat ****************************/

class CohBurgersMat : public CohFrictMat
{
	public :
		virtual ~CohBurgersMat () {};

/// Serialization
	YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohBurgersMat,CohFrictMat,"",
		((Real,kmn,NaN,,"Normal elastic stiffness for Maxwell."))
		((Real,cmn,NaN,,"Normal viscous constant for Maxwell."))
		((Real,kms,NaN,,"Shear elastic stiffness for Maxwell."))
		((Real,cms,NaN,,"Shear viscous constant for Maxwell."))
		((Real,kkn,NaN,,"Normal elastic stiffness for Kelvin."))
		((Real,ckn,NaN,,"Normal viscous constant for Kelvin."))
		((Real,kks,NaN,,"Shear elastic stiffness for Kelvin."))
		((Real,cks,NaN,,"Shear viscous constant for Kelvin.")),
		createIndex();
		);
/// Indexable
	REGISTER_CLASS_INDEX(CohBurgersMat,CohFrictMat);
};

REGISTER_SERIALIZABLE(CohBurgersMat);


//********************** CohBurgersPhys ****************************/

class CohBurgersPhys : public CohFrictPhys
{
	public :
		virtual ~CohBurgersPhys() {};
		void SetBreakingState() {cohesionBroken = true; normalAdhesion = 0; shearAdhesion = 0;};

	YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohBurgersPhys,CohFrictPhys,"",
		((Real,kmn,NaN,,"Stiffness of Maxwell's spring(normal)"))
		((Real,kkn,NaN,,"Stiffness of Kelvin's spring(normal)"))
		((Real,cmn,NaN,,"Viscosity of Maxwell's dashpot(normal)"))
	    	((Real,ckn,NaN,,"Viscosity of Kelvin's dashpot(normal)"))
		((Real,kms,NaN,,"Stiffness of Maxwell's spring(shear)"))
		((Real,kks,NaN,,"Stiffness of Kelvin's spring(shear)"))
		((Real,cms,NaN,,"Viscosity of Maxwell's dashpot(shear)"))
	    	((Real,cks,NaN,,"Viscosity of Kelvin's dashpot(shear)"))
		((Real,MaxDef,NaN,,"maximum permitted deformation"))
		((Vector3r,ukn,Vector3r(0,0,0),,"Normal displacement"))
		((Vector3r,uks,Vector3r(0,0,0),,"Shear displacement")),
		createIndex();
	)
	REGISTER_CLASS_INDEX(CohBurgersPhys,CohFrictPhys);
};

REGISTER_SERIALIZABLE(CohBurgersPhys);


//********************** Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys ****************************/

class Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys: public IPhysFunctor
 {
	public :
		virtual void go(		const shared_ptr<Material>& b1,
						const shared_ptr<Material>& b2,
						const shared_ptr<Interaction>& interaction);
					
		int cohesionDefinitionIteration;				
	YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys,IPhysFunctor,
	"Convert 2 instances of :yref:`CohBurgersMat` to :yref:`CohBurgersPhys` using the rule of consecutive connection.",
	((bool,setCohesionNow,false,,"If true, assign cohesion to all existing contacts"))
	((bool,setCohesionOnNewContacts,false,,"If true, assign cohesion at all new contacts."))
	,
	cohesionDefinitionIteration = -1;
	);
	FUNCTOR2D(CohBurgersMat,CohBurgersMat);
};
REGISTER_SERIALIZABLE(Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys);





class Ip2_CohFrictMat_CohBurgersMat_CohBurgersPhys: public IPhysFunctor
{
	public :
		virtual void go(	const shared_ptr<Material>& b1,
					const shared_ptr<Material>& b2,
					const shared_ptr<Interaction>& interaction);
		int cohesionDefinitionIteration;
	YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_CohFrictMat_CohBurgersMat_CohBurgersPhys,IPhysFunctor,
	"converts an interaction of FrictMat and CohBurgersMat to CohBurgersPhys",
	((bool,setCohesionNow,false,,"If true, assign cohesion to all existing contacts"))
	((bool,setCohesionOnNewContacts,false,,"If true, assign cohesion at all new contacts."))
	,
	cohesionDefinitionIteration = -1;
	);
	FUNCTOR2D(CohFrictMat,CohBurgersMat);
};

REGISTER_SERIALIZABLE(Ip2_CohFrictMat_CohBurgersMat_CohBurgersPhys);




class Ip2_FrictMat_CohBurgersMat_CohBurgersPhys: public IPhysFunctor
{
	public :
		virtual void go(const shared_ptr<Material>& b1,
					const shared_ptr<Material>& b2,
					const shared_ptr<Interaction>& interaction);
		int cohesionDefinitionIteration;
	YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_FrictMat_CohBurgersMat_CohBurgersPhys,IPhysFunctor,
	"converts an interaction of FrictMat and CohBurgersMat to CohBurgersPhys",
	((bool,setCohesionNow,false,,"If true, assign cohesion to all existing contacts"))
	((bool,setCohesionOnNewContacts,false,,"If true, assign cohesion at all new contacts."))
	,
	cohesionDefinitionIteration = -1;
	);
	FUNCTOR2D(FrictMat,CohBurgersMat);
};

REGISTER_SERIALIZABLE(Ip2_FrictMat_CohBurgersMat_CohBurgersPhys);



//******************** Law2_ScGeom6D_CohBurgersPhys_CohesiveBurgers *************************/


class Law2_ScGeom6D_CohBurgersPhys_CohesiveBurgers: public LawFunctor {
	public:
		virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction*);
			
	FUNCTOR2D(ScGeom6D,CohBurgersPhys);
	
	YADE_CLASS_BASE_DOC(Law2_ScGeom6D_CohBurgersPhys_CohesiveBurgers, LawFunctor,"");

};
REGISTER_SERIALIZABLE(Law2_ScGeom6D_CohBurgersPhys_CohesiveBurgers);

bool computeForceTorqueCohBurgers(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r & torque2);