yade-users team mailing list archive

[behzad <question293669@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi folks,

I've got a problem in shear force calculation in a cohesive contact.

That's the case:

Imagine two spheres are in contact like this:

O.bodies.append(sphere([0.0,0.0,0.0],0.01,fixed=True,material=Mat))
O.bodies.append(sphere([0.0,0.0,2.0e-2],0.01,fixed=False,material=Mat))

so the contact normal is in Z direction.

Now imagine the material is cohFrictMat and the contact is cohesive. If a force in Y direction is applied to the upper sphere;

O.forces.addF(1,(0,-1000,0),permanent=True)

the sphere is going to move to left and the contact normal is now parallel to Y axis. 

That's what happens with CohFrictMat and it's normal. 


However, in a new contact model that I have developed, it seems that there's a problem in calculating the shear force that the above mentioned scenario causes the upper sphere to rotate around the lower sphere by an ever increasing speed and the shear force is always increasing. I copy the force calculation code here and I appreciate if one can have a look to see what's the problem.

Thank you.





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();

  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 (geom.penetrationDepth <0 && phys.fragile && phys.normalAdhesion > 0 && phys.normalForce.norm() > phys.normalAdhesion) {
		phys.cohesionBroken= true;
		return false;	
  }
  
   if (phys.fragile && phys.shearAdhesion > 0 && phys.shearForce.norm() > phys.shearAdhesion) {
		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.0*(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 your team yade-users is
an answer contact for Yade.