yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #04705
[Branch ~yade-dev/yade/trunk] Rev 2267: some updates in CohesiveFrictionalPM. Unfortunately, the law still does not work in macroscale si...
------------------------------------------------------------
revno: 2267
committer: Luc Scholtes <sch50p@fluent-ph>
branch nick: trunk
timestamp: Wed 2010-06-02 16:52:37 +1000
message:
some updates in CohesiveFrictionalPM. Unfortunately, the law still does not work in macroscale simulations...
modified:
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.hpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-05-26 08:18:36 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-06-02 06:52:37 +0000
@@ -11,25 +11,27 @@
CREATE_LOGGER(Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM);
void Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* rootBody){
-
+ const Real& dt = rootBody->dt;
+
ScGeom* geom = static_cast<ScGeom*>(ig.get());
CFpmPhys* phys = static_cast<CFpmPhys*>(ip.get());
- int id1 = contact->getId1();
- int id2 = contact->getId2();
- shared_ptr<BodyContainer>& bodies = rootBody->bodies;
- Body* b1 = ( *bodies ) [id1].get();
- Body* b2 = ( *bodies ) [id2].get();
- const Real& dt = scene->dt;
+ const int &id1 = contact->getId1();
+ const int &id2 = contact->getId2();
+ Body* b1 = Body::byId(id1,rootBody).get();
+ Body* b2 = Body::byId(id2,rootBody).get();
Real displN = geom->penetrationDepth; // NOTE: the sign for penetrationdepth is different from ScGeom and Dem3DofGeom: geom->penetrationDepth>0 when spheres interpenetrate
Real Dtensile=phys->FnMax/phys->kn;
Real Dsoftening = phys->strengthSoftening*Dtensile;
+ /// HERE IS THE PROBLEM!!! -> works well for simple cases (2 spheres or sphere/wall interactions), but not for triaxial and uniaxial test -> timestep issue????
/*to set the equilibrium distance between all cohesive elements when they first meet -> stress-free assembly*/
if ( contact->isFresh(rootBody) ) { phys->initD = displN; phys->normalForce = Vector3r::Zero(); phys->shearForce = Vector3r::Zero();}
-
Real D = displN - phys->initD; // interparticular distance is computed depending on the equilibrium distance (which is 0 for non cohesive interactions)
-
+ // those lines make the law equivalent to CohesiveFrictionalContactLaw and make it works in triaxial and uniaxial??????????????????????
+ //if ( contact->isFresh(rootBody) ) { phys->shearForce = Vector3r::Zero();}
+ //Real D = displN;
+
/* Determination of interaction */
if (D < 0){ //spheres do not touch
if (!phys->isCohesive){ rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return; } // destroy the interaction before calculation
@@ -40,6 +42,11 @@
CFpmState* st2=dynamic_cast<CFpmState*>(b2->state.get());
st1->numBrokenCohesive+=1;
st2->numBrokenCohesive+=1;
+ //// the same thing but from ConcretePM
+ //const shared_ptr<Body>& body1=Body::byId(contact->getId1(),rootBody), body2=Body::byId(contact->getId2(),rootBody); assert(body1); assert(body2);
+ //const shared_ptr<CFpmState>& st1=YADE_PTR_CAST<CFpmState>(body1->state), st2=YADE_PTR_CAST<CFpmState>(body2->state);
+ //{ boost::mutex::scoped_lock lock(st1->updateMutex); st1->numBrokenCohesive+=1; }
+ //{ boost::mutex::scoped_lock lock(st2->updateMutex); st2->numBrokenCohesive+=1; }
// end of update
rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return;
}
@@ -60,27 +67,37 @@
/*ShearForce*/
Vector3r& shearForce = phys->shearForce;
- // using scGeom function updateShear(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting)
+ // using scGeom function rotateAndGetShear
State* st1 = Body::byId(id1,rootBody)->state.get();
State* st2 = Body::byId(id2,rootBody)->state.get();
- Vector3r dus = geom->rotateAndGetShear(shearForce, phys->prevNormal, st1, st2, dt, preventGranularRatcheting);
+ // define shift to handle periodicity
+ Vector3r shiftVel = scene->isPeriodic ? (Vector3r)((scene->cell->velGrad*scene->cell->Hsize)*Vector3r((Real) contact->cellDist[0],(Real) contact->cellDist[1],(Real) contact->cellDist[2])) : Vector3r::Zero();
+ Vector3r dus = geom->rotateAndGetShear(shearForce, phys->prevNormal, st1, st2, dt, shiftVel, preventGranularRatcheting);
+ /// before changes to adapt periodic scene it was like that: Vector3r dus = geom->rotateAndGetShear(shearForce, phys->prevNormal, st1, st2, dt, preventGranularRatcheting);
//Linear elasticity giving "trial" shear force
shearForce -= phys->ks*dus;
// needed for the next timestep
phys->prevNormal = geom->normal;
/* Morh-Coulomb criterion */
- Real maxFs = abs(Fn)*phys->tanFrictionAngle + phys->FsMax;
+ Real maxFs = phys->FsMax + Fn*phys->tanFrictionAngle;
if (shearForce.squaredNorm() > maxFs*maxFs){
shearForce*=maxFs/shearForce.norm(); // to fix the shear force to its yielding value
}
/* Apply forces */
- Vector3r f = phys->normalForce + phys->shearForce;
- // NOTE applyForceAtContactPoint computes torque induced by normal and shear force too
- applyForceAtContactPoint(f, geom->contactPoint, contact->getId2(), b2->state->pos, contact->getId1(), b1->state->pos, rootBody); // NOTE id1 and id2 must suit the orientation of the force
-
+ Vector3r f = phys->normalForce + shearForce;
+ // these lines to adapt to periodic boundary conditions (NOTE applyForceAtContactPoint computes torque induced by normal and shear force too)
+ if (!scene->isPeriodic)
+ applyForceAtContactPoint(f , geom->contactPoint , id2, st2->se3.position, id1, st1->se3.position, rootBody);
+ else { // in scg we do not wrap particles positions, hence "applyForceAtContactPoint" cannot be used when scene is periodic
+ rootBody->forces.addForce(id1,-f);
+ rootBody->forces.addForce(id2,f);
+ rootBody->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(-f));
+ rootBody->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(-f));
+ }
+
/* Moment Rotation Law */
// NOTE this part could probably be computed in ScGeom to avoid copy/paste multiplication !!!
Quaternionr delta( b1->state->ori * phys->initialOrientation1.conjugate() *phys->initialOrientation2 * b2->state->ori.conjugate()); //relative orientation
@@ -97,7 +114,7 @@
Vector3r axis_bending(aa.angle()*aa.axis() - axis_twist); //total rotation minus rotation about normal
Vector3r moment_bending(axis_bending * phys->kr);
Vector3r moment = moment_twist + moment_bending;
-
+
Real MomentMax = phys->maxBend*std::fabs(phys->normalForce.norm());
Real scalarMoment = moment.norm();
@@ -116,7 +133,7 @@
rootBody->forces.addTorque(id1,-moment);
rootBody->forces.addTorque(id2, moment);
-
+
}
CREATE_LOGGER(Ip2_CFpmMat_CFpmMat_CFpmPhys);
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.hpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.hpp 2010-04-25 13:18:11 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.hpp 2010-06-02 06:52:37 +0000
@@ -28,7 +28,7 @@
/** This class holds information associated with each body state*/
class CFpmState: public State {
YADE_CLASS_BASE_DOC_ATTRS(CFpmState,State,"CFpm state information about each body.\n\nNone of that is used for computation (at least not now), only for post-processing.",
- ((int,numBrokenCohesive,0,"Number of (cohesive) contacts that damaged completely"))
+ ((int,numBrokenCohesive,0,"Number of broken cohesive links. [-]"))
);
};
REGISTER_SERIALIZABLE(CFpmState);
@@ -106,7 +106,7 @@
FUNCTOR2D(ScGeom,CFpmPhys);
YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM,LawFunctor,"Constitutive law for the CFpm model.",
- ((bool,preventGranularRatcheting,false,"If true rotations are computed such as granular ratcheting is prevented. See article [Alonso2004]_, pg. 3-10 -- and a lot more papers from the same authors)."))
+ ((bool,preventGranularRatcheting,true,"If true rotations are computed such as granular ratcheting is prevented. See article [Alonso2004]_, pg. 3-10 -- and a lot more papers from the same authors)."))
);
DECLARE_LOGGER;
};
