Merge authors: booncw (booncw) ------------------------------------------------------------ revno: 1972 [merge] committer: CWBoon <[email protected]> branch nick: trunk timestamp: Sat 2010-01-16 15:05:03 +0000 message: Cohesionless Moment Rotation LLaw as implemented by Plassiard et al (2009). Verified added: pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.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.
=== added file 'pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp' --- pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp 1970-01-01 00:00:00 +0000 +++ pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp 2010-01-16 15:03:01 +0000 @@ -0,0 +1,277 @@ +/* cwb...@2010 [email protected] */ + + +#include"CohesionlessMomentRotation.hpp" +#include<yade/core/Scene.hpp> +#include<yade/pkg-dem/ScGeom.hpp> +#include<yade/core/Omega.hpp> + +/* Contact law have been verified with Plassiard et al. (2009) : A spherical discrete element model: calibration procedure and incremental response */ + + +YADE_PLUGIN((Law2_SCG_MomentPhys_CohesionlessMomentRotation)(Ip2_BMP_BMP_MomentPhys)(MomentPhys)); + + + +/********************** Law2_Dem3DofGeom_RockPMPhys_Rpm ****************************/ +CREATE_LOGGER(Law2_SCG_MomentPhys_CohesionlessMomentRotation); + +void Law2_SCG_MomentPhys_CohesionlessMomentRotation::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* rootBody){ + + ScGeom* geom = static_cast<ScGeom*>(ig.get()); //InteractionGeometry + MomentPhys* phys = static_cast<MomentPhys*>(ip.get()); //InteractionPhysics + int id1 = contact->getId1(); //Id of body1 + int id2 = contact->getId2(); //Id of body2 + shared_ptr<BodyContainer>& bodies = rootBody->bodies; + Body* b1 = ( *bodies ) [id1].get(); + Body* b2 = ( *bodies ) [id2].get(); + Real dt = Omega::instance().getTimeStep(); //Get TimeStep + + + /*NormalForce */ + Real displN=geom->penetrationDepth; //get penetrationDepth from SCG + if (displN<0){rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return;} //! NOTE: the sign for penetrationdepth is different from SCG and Dem3DofGeom + Real Fn =phys->kn*displN; //scalar force + phys->normalForce= Fn*geom->normal; //vector force NOTE: normal is position2-position1. Therefore it is directed from particle1 to particle2 + + + Vector3r axis; + Real angle; + + + /////////////////////////////////////////////////////////////////////// SHEAR FORCE /////////////////////////////////////////////// + /*ShearForce*/ + Vector3r& shearForce = phys->shearForce; //alias for phys->shearForce + + + /// Here is the code with approximated rotations /// + axis = phys->prevNormal.Cross(geom->normal); //axis of rotation between previous and current + shearForce -= shearForce.Cross(axis); //minus shearForce that is perpendicular (outward) to contact (not tangent) + angle = dt*0.5*geom->normal.Dot(Body::byId(id1)->state->angVel+Body::byId(id2)->state->angVel); + axis = angle*geom->normal; + shearForce -= shearForce.Cross(axis); //minus shearForce that is perpendicular outward to contact (not tangent) + + /// Here is the code with exact rotations /// + //Quaternionr q; + //1st imagine we have to rotate shearForce ON the x-y plane as a result of contact point rotation + //axis = phys->prevNormal.Cross(geom->normal); // axis of rotation of contact point + //angle = Mathr::ACos(geom->normal.Dot(phys->prevNormal)); //angle of rotation of contact point + //q.FromAxisAngle(axis,angle); //quaternion for rotation of contact point + // shearForce = q*shearForce; //rotate shearForce from previous contact orientation to current contact orientation + + //Then imagine we have to find a direction by rotating THROUGH the x-y plane, ABOUT the contact normal (i.e. axis). The direction should be along the resultant angular velocity. + //angle = dt*0.5*geom->normal.Dot(Body::byId(id1)->state->angVel + Body::byId(id1)->state->angVel); //angle of rotation about the contact normal + //axis = geom->normal; + //q.FromAxisAngle(axis,angle); + //shearForce = q*shearForce; + /// /// + + + Vector3r x = geom->contactPoint; + Vector3r c1x = (x - b1->state->pos); //distance of contact point from centre1 + Vector3r c2x = (x - b2 ->state->pos); //distance of contact ponit from centre2 + /// The following definition of c1x and c2x is to avoid "granular ratcheting" + /// (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN, + /// "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena", + /// ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors) + + + Vector3r _c1x_(0,0,0); + Vector3r _c2x_(0,0,0); + + if (preventGranularRatcheting){ + _c1x_ = geom->radius1*geom->normal; //vector radius1 + _c2x_ = -geom->radius2*geom->normal; //vector radius2 + }else{ + _c1x_ = c1x; + _c2x_ = c2x; + } + + Vector3r relativeVelocity = (b2->state->vel + b2->state->angVel.Cross(_c2x_)) - (b1->state->vel + b1->state->angVel.Cross(_c1x_)); + Vector3r shearVelocity = relativeVelocity - geom->normal.Dot(relativeVelocity)*geom->normal; //First find the scalar value of normal velocity. Then change it to a vector. Then minus. + Vector3r shearDisplacement = shearVelocity*dt; //in Dem3DofGeom displacementT() is cumulative + + shearForce -= phys->ks*shearDisplacement; //shearForce is cumulative here because shearDisplacement is not cumulative. See CundallStrack.hpp, it is different! + + + //////////////////////////////////////////////////////////// SHEAR FORCE END /////////////////////////////////////////////////////////////////////////////// + + + /* Coulomb's Law */ + Real maxFs = 0; + Real Fs = shearForce.Length(); //scalar + maxFs = std::max((Real) 0, Fn*phys->tanFrictionAngle); //Fn is scalar + if ( Fs > maxFs ) + { + shearForce *= maxFs/Fs; + //if ( Fn<0 ) phys->normalForce = Vector3r::ZERO; + } + + phys->shearForce = shearForce; //for recording purposes + + /* Apply forces */ + Vector3r f = phys->normalForce + shearForce; + rootBody->forces.addForce (id1,-f); + rootBody->forces.addForce (id2, f); + rootBody->forces.addTorque(id1,-(c1x).Cross(f)); //about axis perpendicular to normal and force + rootBody->forces.addTorque(id2, c2x.Cross(f)); + + /* Moment Rotation Law */ + Quaternionr delta( b1->state->ori * phys->initialOrientation1.Conjugate() *phys->initialOrientation2 * b2->state->ori.Conjugate()); //relative orientation + Vector3r axisM; // axis of rotation - this is the Moment direction UNIT vector. + Real angleM; // angle represents the power of resistant ELASTIC moment + delta.ToAxisAngle(axisM,angleM); + if(angleM > Mathr::PI) angleM -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi + + phys->cumulativeRotation = angleM; + //Real elasticMoment = phys->kr * std::abs(angleM); // positive value (*) + //Vector3r moment = axisM * elasticMoment * (angleM<0.0?-1.0:1.0); // restore sign. (*) + + + //Find angle*axis. That's all. But first find angle about contact normal. Result is scalar. Axis is contact normal. + Real angle_twist(angleM * axisM.Dot(geom->normal) ); //rotation about normal + Vector3r axis_twist(angle_twist * geom->normal); + Vector3r moment_twist(axis_twist * phys->kr); + + Vector3r axis_bending(angleM*axisM - axis_twist); //total rotation minus rotation about normal + Vector3r moment_bending(axis_bending * phys->kr); + + Vector3r moment = moment_twist + moment_bending; + + Real radiusm = 0.5*(geom->radius1 + geom->radius2); //radius is used here + Real MomentMax = phys->Eta * radiusm* std::fabs(Fn); + Real scalarMoment = moment.Length(); + + + if(scalarMoment > MomentMax) /*Plastic moment */ + { + Real ratio=0; + ratio *= MomentMax/scalarMoment; + moment *= ratio; + moment_twist *= ratio; + moment_bending *= ratio; + } + + phys->moment_twist = moment_twist; + phys->moment_bending = moment_bending; + + rootBody->forces.addTorque(id1,-moment); + rootBody->forces.addTorque(id2, moment); + /// Moment law END /// + + phys->prevNormal = geom->normal; + +} + + + +CREATE_LOGGER(Ip2_MomentMat_MomentMat_MomentPhys); + +void Ip2_MomentMat_MomentMat_MomentPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction){ + + if(interaction->interactionPhysics) return; + + ScGeom* scg=dynamic_cast<ScGeom*>(interaction->interactionGeometry.get()); + + assert(scg); + + const shared_ptr<MomentMat>& sdec1 = YADE_PTR_CAST<MomentMat>(b1); + const shared_ptr<MomentMat>& sdec2 = YADE_PTR_CAST<MomentMat>(b2); + + shared_ptr<MomentPhys> contactPhysics(new MomentPhys()); + //interaction->interactionPhysics = shared_ptr<MomentPhys>(new MomentPhys()); + //const shared_ptr<MomentPhys>& contactPhysics = YADE_PTR_CAST<MomentPhys>(interaction->interactionPhysics); + + /* From interaction physics */ + Real Ea = sdec1->young; + Real Eb = sdec2->young; + Real Va = sdec1->poisson; + Real Vb = sdec2->poisson; + Real fa = sdec1->frictionAngle; + Real fb = sdec2->frictionAngle; + Real etaA = sdec1->eta; + Real etaB = sdec2->eta; + + /* From interaction geometry */ + Real Da= scg->radius1; + Real Db= scg->radius2; + + /* calculate stiffness */ + Real Kn=0; + Real Ks=0; + Real Kr=0; + + Real radiusm = (Da + Db)/2; + //Kr = Beta*Ks*radiusm*radiusm; + + if(userInputStiffness == true && useAlphaBeta == true){ + Kn = 0.5*Knormal*radiusm; + Ks = Alpha*Knormal; + Ks = 0.5*Ks*radiusm; + Kr = Beta*Ks*radiusm*radiusm; + } + else if (userInputStiffness == true && useAlphaBeta == false){//if userInputStiffness = true + Kn = Knormal; + Ks = Kshear; + Kr = Krotate; + } + else{ + Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses + Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);//harmonic average of two stiffnesses with ks=V*kn for each sphere + //Kr = Da*Db*Ks*2.0; // just like "2.0" above - it's an arbitrary parameter + Kr=0; + } + + /* Pass values calculated from above to CSPhys */ + contactPhysics->kn = Kn; + contactPhysics->ks = Ks; + contactPhysics->kr = Kr; + contactPhysics->Eta = std::min(etaA,etaB); + contactPhysics->frictionAngle = std::min(fa,fb); + contactPhysics->tanFrictionAngle = std::tan(contactPhysics->frictionAngle); + contactPhysics->initialOrientation1 = Body::byId(interaction->getId1())->state->ori; + contactPhysics->initialOrientation2 = Body::byId(interaction->getId2())->state->ori; + contactPhysics->prevNormal = scg->normal; //This is also done in the Contact Law. It is not redundant because this class is only called ONCE! + + interaction->interactionPhysics = contactPhysics; +} + + + +/* Moment Phys */ +MomentPhys::MomentPhys() +{ + createIndex(); + frictionAngle = 0; + tanFrictionAngle = 0; + Eta = 0; + prevNormal = Vector3r(0,0,0); + + moment_twist = Vector3r(0,0,0); + moment_bending = Vector3r(0,0,0); + + shear = Vector3r(0,0,0); //currently not used + cumulativeRotation =0; +// assign neutral value + initialOrientation1 = Quaternionr(1.0,0.0,0.0,0.0); + initialOrientation2 = Quaternionr(1.0,0.0,0.0,0.0); + kr = 0; +} + +MomentPhys::~MomentPhys(){} + +/* Ip2_BMP_BMP_MomentPhys */ +Ip2_MomentMat_MomentMat_MomentPhys::Ip2_MomentMat_MomentMat_MomentPhys() +{ + userInputStiffness = false; + useAlphaBeta = false; + Alpha = 0; + Beta = 0; + Knormal = 0; + Krotate = 0; + Kshear= 0; + +} +Ip2_MomentMat_MomentMat_MomentPhys::~Ip2_MomentMat_MomentMat_MomentPhys(){}; + + === added file 'pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.hpp' --- pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.hpp 1970-01-01 00:00:00 +0000 +++ pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.hpp 2010-01-16 15:03:01 +0000 @@ -0,0 +1,70 @@ +#pragma once +#include<yade/pkg-common/ElastMat.hpp> +#include<yade/pkg-common/InteractionPhysicsFunctor.hpp> +#include<yade/pkg-common/NormShearPhys.hpp> +#include<yade/pkg-common/LawFunctor.hpp> +#include<yade/pkg-dem/ScGeom.hpp> +#include <set> +#include <boost/tuple/tuple.hpp> + + +class Law2_SCG_MomentPhys_CohesionlessMomentRotation: public LawFunctor{ + public: + virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody); + bool preventGranularRatcheting; + Law2_SCG_MomentPhys_CohesionlessMomentRotation():preventGranularRatcheting(false){}; + FUNCTOR2D(ScGeom,MomentPhys); + REGISTER_CLASS_AND_BASE(Law2_SCG_MomentPhys_CohesionlessMomentRotation,LawFunctor); + REGISTER_ATTRIBUTES(LawFunctor,/*nothing here*/); + DECLARE_LOGGER; +}; +REGISTER_SERIALIZABLE(Law2_SCG_MomentPhys_CohesionlessMomentRotation); + + +class Ip2_MomentMat_MomentMat_MomentPhys: public InteractionPhysicsFunctor{ + public: + bool userInputStiffness, useAlphaBeta; //for users to choose whether to input stiffness directly or use ratios to calculate Ks/Kn + Real Knormal, Kshear, Krotate; //allows user to input stiffness properties from triaxial test. These will be passed to MomentPhys or NormShearPhys + Real Alpha, Beta; //Alpha is a ratio of Ks/Kn, Beta is a ratio to calculate Kr + Ip2_MomentMat_MomentMat_MomentPhys(); + virtual ~Ip2_MomentMat_MomentMat_MomentPhys(); + virtual void go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction); + REGISTER_ATTRIBUTES(InteractionPhysicsFunctor,(userInputStiffness)(useAlphaBeta)(Knormal)(Kshear)(Krotate)(Alpha)(Beta)); + FUNCTOR2D(MomentMat,MomentMat); + REGISTER_CLASS_AND_BASE(Ip2_MomentMat_MomentMat_MomentPhys,InteractionPhysicsFunctor); + DECLARE_LOGGER; +}; +REGISTER_SERIALIZABLE(Ip2_MomentMat_MomentMat_MomentPhys); + + +class MomentPhys: public NormShearPhys { + private: + public: + Real frictionAngle, tanFrictionAngle, Eta; + Vector3r prevNormal; + Vector3r shear; + Real cumulativeRotation; + Quaternionr initialOrientation1,initialOrientation2; + Real kr; // rolling stiffness + Vector3r moment_twist,moment_bending; + MomentPhys(); + + virtual ~MomentPhys(); + + REGISTER_ATTRIBUTES(NormShearPhys,(tanFrictionAngle) (frictionAngle) (prevNormal) (initialOrientation1) (initialOrientation2) (kr) (Eta)(moment_twist) (moment_bending) (shear)(cumulativeRotation)); + REGISTER_CLASS_AND_BASE(MomentPhys,NormShearPhys); + REGISTER_CLASS_INDEX(MomentPhys,NormShearPhys); +}; +REGISTER_SERIALIZABLE(MomentPhys); + +/** This class holds information associated with each body */ +class MomentMat: public FrictMat { + public: + Real eta; //It has to be stored in this class and not by ContactLaw, because users may want to change its values before/after isotropic compaction. + MomentMat(): eta(0) {createIndex();}; + REGISTER_ATTRIBUTES(FrictMat, (eta)); + REGISTER_CLASS_AND_BASE(MomentMat,FrictMat); + REGISTER_CLASS_INDEX(MomentMat,FrictMat); +}; +REGISTER_SERIALIZABLE(MomentMat); +
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
