yade-users team mailing list archive

[Jan Stránský <question293669@xxxxxxxxxxxxxxxxxxxxx>]

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

Jan Stránský proposed the following answer:
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ý <
question293669@xxxxxxxxxxxxxxxxxxxxx>:

> 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 <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.
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~yade-users
> > Post to     : yade-users@xxxxxxxxxxxxxxxxxxx
> > 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     : yade-users@xxxxxxxxxxxxxxxxxxx
> 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.