yade-users team mailing list archive

[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.