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ý < [email protected]>: > 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 protected]>: > > > 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 protected] > > 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 protected] > 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 protected] Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
