yade-users team mailing list archive

[behzad <question261724@xxxxxxxxxxxxxxxxxxxxx>]

New question #261724 on Yade:
https://answers.launchpad.net/yade/+question/261724

Hi guys,

I created a new class of material to present viscoelasticity by Burgers model. This class of material has been derived from CohFrictMat and it's called CohBurgersMat. I verified the response of CohBurgersMat-CohBurgersMat interaction by comparing the load-displacement response of analytical equations. 

However, I need to have FrictMat-CohBurgersMat interaction as well. To do so, I created a new Ip2 (Ip2_FrictMat_CohBurgersMat_CohBurgersPhys).

However, the model fails to create the required interaction physics. Can you please have a look on my script? Thanks.



======== CohBurgers.cpp ==============================================

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

#include"CohBurgers.hpp"
#include<core/State.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_FrictMat_CohBurgersMat_CohBurgersPhys)(Law2_ScGeom_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());
  Calculate_FrictMat_CohBurgersMat_CohBurgersPhys(b1, b2, interaction, phys);
  interaction->phys = phys;

  phys->cohesionBroken = false;
  phys->initCohesion   = true;
  cohesionDefinitionIteration = -1;
  setCohesionNow = true;
}

//********************** 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());
  Calculate_CohBurgersMat_CohBurgersMat_CohBurgersPhys(b1, b2, interaction, phys);
  interaction->phys = phys;
  
  phys->cohesionBroken = false;
  phys->initCohesion   = true;
  cohesionDefinitionIteration = -1;
  setCohesionNow = true;
}


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

bool Law2_ScGeom_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=*static_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;
}

/*!Calculate_CohBurgersPhys (1)*/

void Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys::Calculate_CohBurgersMat_CohBurgersMat_CohBurgersPhys(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction, shared_ptr<CohBurgersPhys> phys) {

  if(interaction->phys) return;
  CohBurgersMat* mat1 = static_cast<CohBurgersMat*>(b1.get());
  CohBurgersMat* mat2 = static_cast<CohBurgersMat*>(b2.get());
  Real mass1 = 1.0;
  Real mass2 = 1.0;
	
  mass1=Body::byId(interaction->getId1())->state->mass;
  mass2=Body::byId(interaction->getId2())->state->mass;

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

  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);
  phys->normalAdhesion= (mat1->normalCohesion);

  interaction->phys = shared_ptr<CohBurgersPhys>(phys);
}


/*!Calculate_CohBurgersPhys (2)*/

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

  if(interaction->phys) return;
  const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1);
  const shared_ptr<CohBurgersMat>& mat2 = YADE_PTR_CAST<CohBurgersMat>(b2);

  Real mass1 = 1.0;
  Real mass2 = 1.0;
	
  mass1=Body::byId(interaction->getId1())->state->mass;
  mass2=Body::byId(interaction->getId2())->state->mass;

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

  phys->kmn = contactParameterCalculation(kn1,kmn2);
  phys->kms = contactParameterCalculation(ks1,kms2);
  phys->cmn = cmn2;
  phys->cms = cms2;

  phys->kkn = kkn2;
  phys->kks = kks2;
  phys->ckn = ckn2;
  phys->cks = cks2;
  phys->normalAdhesion= ((mat2->normalCohesion));

  interaction->phys = shared_ptr<CohBurgersPhys>(phys);
}


/*! computeForceTorqueBurgers */
 

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

  const ScGeom& geom=*static_cast<ScGeom*>(_geom.get());
  CohBurgersPhys& phys=*static_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 = *static_cast<State*>(bodies[id1]->state.get());
  const State& de2 = *static_cast<State*>(bodies[id2]->state.get());

  const Real& dt = scene->dt;
 
    if (phys.normalForce.norm() > phys.normalAdhesion) {
		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;
    } 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;   
      
  }

}

==========================================================================
 

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.