Yade

shear Force in a cohesive contact

Asked by behzad

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;

  }

}

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said : 		#1

Hi Behzad,
have a look to other Law2::go, you should have something like
geom->rotate(shearForce) and/or geom->rotate(uks) and probably also
geom->precompute there..
cheers
Jan

2016-05-13 16:52 GMT+02:00 behzad <email address hidden>:

> 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.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
Jan Stránský (honzik) said : 		#2

Sorry, precompute is matter of geom itself, so it should not be used in
Law2, but rotate is definitely missing.
Jan

2016-05-13 17:07 GMT+02:00 Jan Stránský <
<email address hidden>>:

> Question #293669 on Yade changed:
> https://answers.launchpad.net/yade/+question/293669
>
> Status: Open => Answered
>
> Jan Stránský proposed the following answer:
> Hi Behzad,
> have a look to other Law2::go, you should have something like
> geom->rotate(shearForce) and/or geom->rotate(uks) and probably also
> geom->precompute there..
> cheers
> Jan
>
>
> 2016-05-13 16:52 GMT+02:00 behzad <email address hidden>:
>
> > 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.
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~yade-users
> > Post to : <email address hidden>
> > Unsubscribe : https://launchpad.net/~yade-users
> > More help : https://help.launchpad.net/ListHelp
> >
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
behzad (behzad-majidi) said : 		#3

I included

phys.shearForce=geom.rotate(phys.shearForce); in the Law2 but it doesn't change anything!

      phys.shearForce=geom.rotate(phys.shearForce);
      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;

Revision history for this message
behzad (behzad-majidi) said : 		#4

one update:

if I remove the following piece of the code, it behaves as expected.

      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;
      }

AFAIK, this piece controls the shear force. Do I need this is this type of contact?
By applying the force as explained in #1, the contact cohesion is not being broken, the force only moves the upper sphere and so contact point is moving.

Revision history for this message
behzad (behzad-majidi) said : 		#5

second update:

the complete normal situation is when I use both:

geom->rotate(shearForce)

and

geom->rotate(normalForce)

Is this normal?

Revision history for this message
Jan Stránský (honzik) said : 		#6

Hi Behzad,

      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;

these values should be available in IGeom class, you don't need to compute
it "manually". Code duplication is always possible source of problems.

anyway, I think now the best way how to find the problem is to compare your
law to some existing law, which works ok (CohFrictMat?).

Another option is some debugging, checking actual values and comparing them
to expected/correct/analytical values)..

cheers
Jan

2016-05-14 2:17 GMT+02:00 behzad <email address hidden>:

> Question #293669 on Yade changed:
> https://answers.launchpad.net/yade/+question/293669
>
> behzad gave more information on the question:
>
> second update:
>
> the complete normal situation is when I use both:
>
> geom->rotate(shearForce)
>
> and
>
> geom->rotate(normalForce)
>
> Is this normal?
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Can you help with this problem?

Provide an answer of your own, or ask behzad for more information if necessary.

To post a message you must log in.