yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06690
[Branch ~yade-dev/yade/trunk] Rev 2646: - Get rid of prevNormal (MindlinPhys inherits again from FrictPhys)
------------------------------------------------------------
revno: 2646
committer: Chiara Modenese <c.modenese@xxxxxxxxx>
branch nick: yade
timestamp: Tue 2011-01-11 19:01:11 +0000
message:
- Get rid of prevNormal (MindlinPhys inherits again from FrictPhys)
- Small changes to HM (I added the code for the moment rotation but using quaternions, work in progress to use the incremental formulation instead of that)
modified:
pkg/dem/HertzMindlin.cpp
pkg/dem/HertzMindlin.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/HertzMindlin.cpp'
--- pkg/dem/HertzMindlin.cpp 2011-01-02 10:12:16 +0000
+++ pkg/dem/HertzMindlin.cpp 2011-01-11 19:01:11 +0000
@@ -1,4 +1,4 @@
-// 2010 © Chiara Modenese <c.modenese@xxxxxxxxx>
+// 2010 © Chiara Modenese <c.modenese@xxxxxxxxx>
#include"HertzMindlin.hpp"
@@ -6,6 +6,7 @@
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
+
YADE_PLUGIN(
(MindlinPhys)
(Ip2_FrictMat_FrictMat_MindlinPhys)
@@ -34,8 +35,9 @@
interaction->phys = mindlinPhys;
FrictMat* mat1 = YADE_CAST<FrictMat*>(b1.get());
FrictMat* mat2 = YADE_CAST<FrictMat*>(b2.get());
-
-
+
+
+
/* from interaction physics */
Real Ea = mat1->young;
Real Eb = mat2->young;
@@ -59,6 +61,7 @@
Real V = (Va+Vb)/2; // average of poisson's ratio
Real E = Ea*Eb/((1.-std::pow(Va,2))*Eb+(1.-std::pow(Vb,2))*Ea); // Young modulus
Real R = Da*Db/(Da+Db); // equivalent radius
+ Real Rmean = (Da+Db)/2.; // mean radius
Real Kno = 4./3.*E*sqrt(R); // coefficient for normal stiffness
Real Kso = 2*sqrt(4*R)*G/(2-V); // coefficient for shear stiffness
Real frictionAngle = std::min(fa,fb);
@@ -66,33 +69,31 @@
Real Adhesion = 4.*Mathr::PI*R*gamma; // calculate adhesion force as predicted by DMT theory
/* pass values calculated from above to MindlinPhys */
- mindlinPhys->tangensOfFrictionAngle = std::tan(frictionAngle);
- mindlinPhys->prevNormal = scg->normal;
+ mindlinPhys->tangensOfFrictionAngle = std::tan(frictionAngle);
mindlinPhys->kno = Kno; // this is just a coeff
mindlinPhys->kso = Kso; // this is just a coeff
mindlinPhys->adhesionForce = Adhesion;
+
+ mindlinPhys->kr = krot;
+ mindlinPhys->maxBendPl = eta*Rmean; // does this make sense? why do we take Rmean?
/* compute viscous coefficients */
if(en && betan) throw std::invalid_argument("Ip2_FrictMat_FrictMat_MindlinPhys: only one of en, betan can be specified.");
if(es && betas) throw std::invalid_argument("Ip2_FrictMat_FrictMat_MindlinPhys: only one of es, betas can be specified.");
- // en specified, compute betan
- if(en){
- Real logE=log((*en)(mat1->id,mat2->id));
- mindlinPhys->betan=-logE/sqrt(pow(Mathr::PI,2)+pow(logE,2));
+ // en or es specified, just compute alpha, otherwise alpha remains 0
+ if(en || es){
+ Real logE = log((*en)(mat1->id,mat2->id));
+ mindlinPhys->alpha = -sqrt(5/6.)*2*logE/sqrt(pow(logE,2)+pow(Mathr::PI,2))*sqrt(2*E*sqrt(R)); // (see Tsuji, 1992)
}
+
// betan specified, use that value directly; otherwise give zero
- else mindlinPhys->betan=betan ? (*betan)(mat1->id,mat2->id) : 0;
-
- // same for es, except that when not given, use the value of betan
- if(es){
- Real logE=log((*es)(mat1->id,mat2->id));
- mindlinPhys->betas=-logE/sqrt(pow(Mathr::PI,2)+pow(logE,2));
+ else{
+ mindlinPhys->betan=betan ? (*betan)(mat1->id,mat2->id) : 0;
+ mindlinPhys->betas=betas ? (*betas)(mat1->id,mat2->id) : mindlinPhys->betan;
}
- else mindlinPhys->betas=betas ? (*betas)(mat1->id,mat2->id) : mindlinPhys->betan;
}
-
/* Function to count the number of adhesive contacts in the simulation at each time step */
Real Law2_ScGeom_MindlinPhys_Mindlin::contactsAdhesive() // It is returning something rather than zero only if includeAdhesion is set to true
{
@@ -111,7 +112,7 @@
Real normEnergy=0;
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
if(!I->isReal()) continue;
- ScGeom* scg = dynamic_cast<ScGeom*>(I->geom.get());
+ ScGeom6D* scg = dynamic_cast<ScGeom6D*>(I->geom.get());
MindlinPhys* phys = dynamic_cast<MindlinPhys*>(I->phys.get());
if (phys) {
if (includeAdhesion) {normEnergy += (std::pow(scg->penetrationDepth,5./2.)*2./5.*phys->kno - phys->adhesionForce*scg->penetrationDepth);}
@@ -121,10 +122,27 @@
return normEnergy;
}
+/* Function to get the adhesion energy of the system */
+Real Law2_ScGeom_MindlinPhys_Mindlin::adhesionEnergy()
+{
+ Real adhesionEnergy=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ ScGeom6D* scg = dynamic_cast<ScGeom6D*>(I->geom.get());
+ MindlinPhys* phys = dynamic_cast<MindlinPhys*>(I->phys.get());
+ if (phys && includeAdhesion) {
+ Real R=scg->radius1*scg->radius2/(scg->radius1+scg->radius2);
+ Real gammapi=phys->adhesionForce/(4.*R);
+ adhesionEnergy += gammapi*pow(phys->radius,2);} // note that contact radius is calculated if we calculate energy components
+ }
+ return adhesionEnergy;
+}
+
+
void Law2_ScGeom_MindlinPhys_MindlinDeresiewitz::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
Body::id_t id1(contact->getId1()), id2(contact->getId2());
ScGeom* geom = static_cast<ScGeom*>(ig.get());
- MindlinPhys* phys=static_cast<MindlinPhys*>(ip.get());
+ MindlinPhys* phys=static_cast<MindlinPhys*>(ip.get());
const Real uN=geom->penetrationDepth;
if (uN<0) {scene->interactions->requestErase(id1,id2); return;}
// normal force
@@ -137,7 +155,7 @@
// contact radius
Real R=geom->radius1*geom->radius2/(geom->radius1+geom->radius2);
phys->radius=pow(Fn*pow(R,3/2.)/phys->kno,1/3.);
-
+
// shear force: transform, but keep the old value for now
geom->rotate(phys->usTotal);
Vector3r usOld=phys->usTotal;
@@ -153,7 +171,7 @@
if(Fs.squaredNorm()>maxFs2) Fs*=sqrt(maxFs2)/Fs.norm();
#endif
// apply forces
- Vector3r f=-phys->normalForce-phys->shearForce;
+ Vector3r f=-phys->normalForce-phys->shearForce;
scene->forces.addForce(id1,f);
scene->forces.addForce(id2,-f);
scene->forces.addTorque(id1,(geom->radius1-.5*geom->penetrationDepth)*geom->normal.cross(f));
@@ -163,23 +181,25 @@
void Law2_ScGeom_MindlinPhys_HertzWithLinearShear::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
Body::id_t id1(contact->getId1()), id2(contact->getId2());
ScGeom* geom = static_cast<ScGeom*>(ig.get());
- MindlinPhys* phys=static_cast<MindlinPhys*>(ip.get());
+ MindlinPhys* phys=static_cast<MindlinPhys*>(ip.get());
const Real uN=geom->penetrationDepth;
if (uN<0) {scene->interactions->requestErase(id1,id2); return;}
// normal force
Real Fn=phys->kno*pow(uN,3/2.);
phys->normalForce=Fn*geom->normal;
//phys->kn=3./2.*phys->kno*std::pow(uN,0.5); // update stiffness, not needed
-
+
// shear force
Vector3r& Fs=geom->rotate(phys->shearForce);
Real ks= nonLin>0 ? phys->kso*std::pow(uN,0.5) : phys->kso;
Vector3r shearIncrement;
if(nonLin>1){
- State *de1=Body::byId(id1,scene)->state.get(), *de2=Body::byId(id2,scene)->state.get();
- Vector3r shift2=scene->isPeriodic ? Vector3r(scene->cell->Hsize*contact->cellDist.cast<Real>()) : Vector3r::Zero();
+ State *de1=Body::byId(id1,scene)->state.get(), *de2=Body::byId(id2,scene)->state.get();
Vector3r shiftVel=scene->isPeriodic ? Vector3r(scene->cell->velGrad*scene->cell->Hsize*contact->cellDist.cast<Real>()) : Vector3r::Zero();
- Vector3r incidentV = geom->getIncidentVel(de1, de2, scene->dt, shift2, shiftVel, /*preventGranularRatcheting*/ nonLin>2 );
+ Vector3r shift2 = scene->isPeriodic ? Vector3r(scene->cell->Hsize*contact->cellDist.cast<Real>()): Vector3r::Zero();
+
+
+ Vector3r incidentV = geom->getIncidentVel(de1, de2, scene->dt, shift2, shiftVel, /*preventGranularRatcheting*/ nonLin>2 );
Vector3r incidentVn = geom->normal.dot(incidentV)*geom->normal; // contact normal velocity
Vector3r incidentVs = incidentV-incidentVn; // contact shear velocity
shearIncrement=incidentVs*scene->dt;
@@ -190,7 +210,7 @@
if(Fs.squaredNorm()>maxFs2) Fs*=sqrt(maxFs2)/Fs.norm();
// apply forces
- Vector3r f=-phys->normalForce-phys->shearForce; /* should be a reference returned by geom->rotate */ assert(phys->shearForce==Fs);
+ Vector3r f=-phys->normalForce-phys->shearForce; /* should be a reference returned by geom->rotate */ assert(phys->shearForce==Fs);
scene->forces.addForce(id1,f);
scene->forces.addForce(id2,-f);
scene->forces.addTorque(id1,(geom->radius1-.5*geom->penetrationDepth)*geom->normal.cross(f));
@@ -203,20 +223,21 @@
void Law2_ScGeom_MindlinPhys_Mindlin::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
const Real& dt = scene->dt; // get time step
-
+
Body::id_t id1 = contact->getId1(); // get id body 1
Body::id_t id2 = contact->getId2(); // get id body 2
State* de1 = Body::byId(id1,scene)->state.get();
- State* de2 = Body::byId(id2,scene)->state.get();
-
- ScGeom* scg = static_cast<ScGeom*>(ig.get());
- MindlinPhys* phys = static_cast<MindlinPhys*>(ip.get());
-
- const shared_ptr<Body>& b1=Body::byId(id1,scene); // not sure here (I need it because I want to know if the body isDynamic or not)
- const shared_ptr<Body>& b2=Body::byId(id2,scene); // not sure here
-
- bool useDamping=(phys->betan!=0. || phys->betas!=0.);
+ State* de2 = Body::byId(id2,scene)->state.get();
+
+ ScGeom6D* scg = static_cast<ScGeom6D*>(ig.get());
+ MindlinPhys* phys = static_cast<MindlinPhys*>(ip.get());
+
+ const shared_ptr<Body>& b1=Body::byId(id1,scene);
+ const shared_ptr<Body>& b2=Body::byId(id2,scene);
+
+ bool useDamping=(phys->betan!=0. || phys->betas!=0. || phys->alpha!=0.);
+ if(phys->alpha==0. && !LinDamp) throw std::invalid_argument("Law2_ScGeom_MindlinPhys_Mindlin: viscous contact damping for non-linear elastic law is specified, specify en and es in Ip2.");
// tangential and normal stiffness coefficients, recomputed from betan,betas at every step
Real cn=0, cs=0;
@@ -224,10 +245,10 @@
/****************/
/* NORMAL FORCE */
/****************/
-
- Real uN = scg->penetrationDepth; // get overlapping
+
+ Real uN = scg->penetrationDepth; // get overlapping
if (uN<0) {scene->interactions->requestErase(contact->getId1(),contact->getId2()); return;}
- /* Hertz-Mindlin's formulation (PFC)
+ /* Hertz-Mindlin's formulation (PFC)
Note that the normal stiffness here is a secant value (so as it is cannot be used in the GSTS)
In the first place we get the normal force and then we store kn to be passed to the GSTS */
Real Fn = phys->kno*std::pow(uN,1.5); // normal Force (scalar)
@@ -239,27 +260,28 @@
if (calcEnergy){
Real R=scg->radius1*scg->radius2/(scg->radius1+scg->radius2);
- phys->radius=pow((Fn+(includeAdhesion?phys->adhesionForce:0.))*pow(R,3/2.)/phys->kno,1/3.); // not used anywhere
+ phys->radius=pow((Fn+(includeAdhesion?phys->adhesionForce:0.))*pow(R,3/2.)/phys->kno,1/3.); // attribute not used anywhere, we do not need it
}
/*******************************/
/* TANGENTIAL NORMAL STIFFNESS */
/*******************************/
-
+
phys->kn = 3./2.*phys->kno*std::pow(uN,0.5); // here we store the value of kn to compute the time step
-
+
/******************************/
/* TANGENTIAL SHEAR STIFFNESS */
/******************************/
-
+
phys->ks = phys->kso*std::pow(uN,0.5); // get tangential stiffness (this is a tangent value, so we can pass it to the GSTS)
/************************/
/* DAMPING COEFFICIENTS */
/************************/
-
+
// Inclusion of local damping if requested
- if (useDamping){
+ // viscous damping is defined for both linear and non-linear elastic case
+ if (useDamping && LinDamp){
Real mbar = (!b1->isDynamic() && b2->isDynamic()) ? de2->mass : ((!b2->isDynamic() && b1->isDynamic()) ? de1->mass : (de1->mass*de2->mass / (de1->mass + de2->mass))); // get equivalent mass if both bodies are dynamic, if not set it equal to the one of the dynamic body
//Real mbar = de1->mass*de2->mass / (de1->mass + de2->mass); // equivalent mass
Real Cn_crit = 2.*sqrt(mbar*phys->kn); // Critical damping coefficient (normal direction)
@@ -269,33 +291,36 @@
cs = Cs_crit*phys->betas; // Damping tangential coefficient
if(phys->kn<0 || phys->ks<0){ cerr<<"Negative stiffness kn="<<phys->kn<<" ks="<<phys->ks<<" for ##"<<b1->getId()<<"+"<<b2->getId()<<", step "<<scene->iter<<endl; }
}
+ else if (useDamping){ // (see Tsuji, 1992)
+ Real mbar = (!b1->isDynamic() && b2->isDynamic()) ? de2->mass : ((!b2->isDynamic() && b1->isDynamic()) ? de1->mass : (de1->mass*de2->mass / (de1->mass + de2->mass))); // get equivalent mass if both bodies are dynamic, if not set it equal to the one of the dynamic body
+ cn = phys->alpha*sqrt(mbar)*pow(uN,0.25); // normal viscous coefficient
+ cs = cn; // same value for shear viscous coefficient
+ }
/***************/
/* SHEAR FORCE */
/***************/
-
+
Vector3r& shearElastic = phys->shearElastic; // reference for shearElastic force
- // Define shift to handle periodicity
- Vector3r shift2 = scene->isPeriodic ? Vector3r( scene->cell->Hsize *contact->cellDist.cast<Real>()) : Vector3r::Zero();
- Vector3r shiftVel = scene->isPeriodic ? Vector3r((scene->cell->velGrad*scene->cell->Hsize)*contact->cellDist.cast<Real>()) : Vector3r::Zero();
+ // Define shifts to handle periodicity
+ const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero();
+ const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(contact->cellDist): Vector3r::Zero();
// 1. Rotate shear force
shearElastic = scg->rotate(shearElastic);
Vector3r prev_FsElastic = shearElastic; // save shear force at previous time step
// 2. Get incident velocity, get shear and normal components
- // FIXME: Concerning the possibility to avoid granular ratcheting, it is not clear how this works in the case of HM. See the thread http://www.mail-archive.com/yade-users@xxxxxxxxxxxxxxxxxxx/msg01947.html
- Vector3r incidentV = scg->getIncidentVel(de1, de2, dt, shift2, shiftVel, preventGranularRatcheting);
+ Vector3r incidentV = scg->getIncidentVel(de1, de2, dt, shift2, shiftVel, preventGranularRatcheting);
Vector3r incidentVn = scg->normal.dot(incidentV)*scg->normal; // contact normal velocity
Vector3r incidentVs = incidentV - incidentVn; // contact shear velocity
// 3. Get shear force (incrementally)
shearElastic = shearElastic - phys->ks*(incidentVs*dt);
-
/**************************************/
/* VISCOUS DAMPING (Normal direction) */
/**************************************/
-
+
// normal force must be updated here before we apply the Mohr-Coulomb criterion
- if (useDamping){ // get normal viscous component
+ if (useDamping){ // get normal viscous component
phys->normalViscous = cn*incidentVn;
// add normal viscous component if damping is included
phys->normalForce -= phys->normalViscous;
@@ -305,7 +330,7 @@
/*************************************/
/* SHEAR DISPLACEMENT (elastic only) */
/*************************************/
-
+
Vector3r& us_elastic = phys->usElastic;
us_elastic = scg->rotate(us_elastic); // rotate vector
Vector3r prevUs_el = us_elastic; // store previous elastic shear displacement (already rotated)
@@ -314,19 +339,18 @@
/****************************************/
/* SHEAR DISPLACEMENT (elastic+plastic) */
/****************************************/
-
+
Vector3r& us_total = phys->usTotal;
us_total = scg->rotate(us_total); // rotate vector
Vector3r prevUs_tot = us_total; // store previous total shear displacement (already rotated)
us_total -= incidentVs*dt; // add shear increment NOTE: this vector is not passed into the failure criterion, hence it holds also the plastic part of the shear displacement
-
bool noShearDamp = false; // bool to decide whether we need to account for shear damping dissipation or not
-
+
/********************/
/* MOHR-COULOMB law */
/********************/
-
+
phys->shearViscous=Vector3r::Zero(); // reset so that during sliding, the previous values is not there
if (!includeAdhesion) {
Real maxFs = Fn*phys->tangensOfFrictionAngle;
@@ -357,7 +381,7 @@
/************************/
/* SHEAR ELASTIC ENERGY */
/************************/
-
+
// NOTE: shear elastic energy calculation must come after the MC criterion, otherwise displacements and forces are not updated
if (calcEnergy) {
shearEnergy += (us_elastic-prevUs_el).dot((shearElastic+prev_FsElastic)/2.); // NOTE: no additional energy if we perform sliding since us_elastic and prevUs_el will hold the same value (in fact us_elastic is only keeping the elastic part). We work out the area of the trapezium.
@@ -366,7 +390,7 @@
/**************************************************/
/* VISCOUS DAMPING (energy term, shear direction) */
/**************************************************/
-
+
if (useDamping){ // get normal viscous component (the shear one is calculated inside Mohr-Coulomb criterion, see above)
if (calcEnergy) {if (!noShearDamp) {shearDampDissip += phys->shearViscous.dot(incidentVs*dt);}} // calc energy dissipation due to viscous linear damping
}
@@ -374,7 +398,7 @@
/****************/
/* APPLY FORCES */
/****************/
-
+
if (!scene->isPeriodic)
applyForceAtContactPoint(-phys->normalForce - phys->shearForce, scg->contactPoint , id1, de1->se3.position, id2, de2->se3.position);
else { // in scg we do not wrap particles positions, hence "applyForceAtContactPoint" cannot be used
@@ -384,61 +408,24 @@
scene->forces.addTorque(id1,(scg->radius1-0.5*scg->penetrationDepth)* scg->normal.cross(force));
scene->forces.addTorque(id2,(scg->radius2-0.5*scg->penetrationDepth)* scg->normal.cross(force));
}
- phys->prevNormal = scg->normal;
-}
-
-
-
-#if 0
-void Law2_ScGeom_MindlinPhys_Mindlin::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact, Scene* ncb){
- const Real& dt = scene->dt; // get time step
-
- int id1 = contact->getId1(); // get id body 1
- int id2 = contact->getId2(); // get id body 2
-
- State* de1 = Body::byId(id1,ncb)->state.get();
- State* de2 = Body::byId(id2,ncb)->state.get();
-
- ScGeom* scg = static_cast<ScGeom*>(ig.get());
- MindlinPhys* phys = static_cast<MindlinPhys*>(ip.get());
-
-
- /*** NORMAL FORCE ****/
- Real uN = scg->penetrationDepth; // get overlapping
- if (uN<0) {ncb->interactions->requestErase(contact->getId1(),contact->getId2()); return;}
- /*** Hertz-Mindlin's formulation (PFC) ***/
- phys->kn = phys->kno*sqrt(uN); // normal stiffness
- Real Fn = phys->kn*uN; // normal Force (scalar)
- phys->normalForce = Fn*scg->normal; // normal Force (vector)
-
-
- /*** SHEAR FORCE ***/
- phys->ks = phys->kso*std::pow(Fn,1/3); // shear stiffness
- Vector3r& trialFs = phys->shearForce;
- 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 shearDisp = scg->rotateAndGetShear(trialFs,phys->prevNormal,de1,de2,dt,shiftVel,preventGranularRatcheting);
- trialFs -= phys->ks*shearDisp;
-
-
- /*** MOHR-COULOMB LAW ***/
- Real maxFs = phys->normalForce.squaredNorm();
- if (trialFs.squaredNorm() > maxFs)
- {Real ratio = sqrt(maxFs)/trialFs.norm(); trialFs *= ratio;}
-
-
- /*** APPLY FORCES ***/
- if (!scene->isPeriodic)
- applyForceAtContactPoint(-phys->normalForce - trialFs, scg->contactPoint , id1, de1->se3.position, id2, de2->se3.position, ncb);
- else { // in scg we do not wrap particles positions, hence "applyForceAtContactPoint" cannot be used
- Vector3r force = -phys->normalForce - trialFs;
- ncb->forces.addForce(id1,force);
- ncb->forces.addForce(id2,-force);
- ncb->forces.addTorque(id1,(scg->radius1-0.5*scg->penetrationDepth)* scg->normal.cross(force));
- ncb->forces.addTorque(id2,(scg->radius2-0.5*scg->penetrationDepth)* scg->normal.cross(force));
+
+ /********************************************/
+ /* MOMENT CONTACT LAW (only bending moment) */
+ /********************************************/
+
+ if (includeMoment){
+ phys->moment_bending = scg->getBending()*phys->kr;
+ Real MomentMax = phys->maxBendPl*phys->normalForce.norm();
+ Real scalarMoment = phys->moment_bending.norm();
+ if(scalarMoment > MomentMax)
+ {
+ Real ratio = MomentMax/scalarMoment; // to fix the moment to its yielding value
+ phys->moment_bending *= ratio;
+ }
+ scene->forces.addTorque(id1,-phys->moment_bending);
+ scene->forces.addTorque(id2, phys->moment_bending);
}
- phys->prevNormal = scg->normal;
}
-#endif
=== modified file 'pkg/dem/HertzMindlin.hpp'
--- pkg/dem/HertzMindlin.hpp 2011-01-02 10:12:16 +0000
+++ pkg/dem/HertzMindlin.hpp 2011-01-11 19:01:11 +0000
@@ -1,9 +1,10 @@
-// 2010 © Chiara Modenese <c.modenese@xxxxxxxxx>
-//
+// 2010 © Chiara Modenese <c.modenese@xxxxxxxxx>
+//
/*
=== HIGH LEVEL OVERVIEW OF Mindlin ===
Mindlin is a set of classes to include the Hertz-Mindlin formulation for the contact stiffnesses.
+The DMT formulation is also considered (for adhesive particles, rigid and small bodies).
*/
@@ -24,37 +25,40 @@
/******************** MindlinPhys *********************************/
-class MindlinPhys: public FrictPhysTransitory{
+class MindlinPhys: public FrictPhys{
public:
virtual ~MindlinPhys();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(MindlinPhys,FrictPhysTransitory,"Representation of an interaction of the Hertz-Mindlin type.",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(MindlinPhys,FrictPhys,"Representation of an interaction of the Hertz-Mindlin type.",
((Real,kno,0.0,,"Constant value in the formulation of the normal stiffness"))
((Real,kso,0.0,,"Constant value in the formulation of the tangential stiffness"))
+ ((Real,kr,0.0,,"Rotational stiffness"))
+ ((Real,maxBendPl,0.0,,"Coefficient to determine the maximum plastic moment to apply at the contact"))
+
((Vector3r,normalViscous,Vector3r::Zero(),,"Normal viscous component"))
((Vector3r,shearViscous,Vector3r::Zero(),,"Shear viscous component"))
((Vector3r,shearElastic,Vector3r::Zero(),,"Total elastic shear force"))
((Vector3r,usElastic,Vector3r::Zero(),,"Total elastic shear displacement (only elastic part)"))
((Vector3r,usTotal,Vector3r::Zero(),,"Total elastic shear displacement (elastic+plastic part)"))
+ ((Vector3r,moment_bending,Vector3r::Zero(),,"Artificial bending moment to provide rolling resistance in order to account for some degree of interlocking between particles"))
+
((Real,radius,NaN,,"Contact radius (only computed with :yref:`Law2_ScGeom_MindlinPhys_Mindlin::calcEnergy`)"))
//((Real,gamma,0.0,"Surface energy parameter [J/m^2] per each unit contact surface, to derive DMT formulation from HM"))
((Real,adhesionForce,0.0,,"Force of adhesion as predicted by DMT"))
((bool,isAdhesive,false,,"bool to identify if the contact is adhesive, that is to say if the contact force is attractive"))
+ // Contact damping coefficients as for linear elastic contact law
((Real,betan,0.0,,"Fraction of the viscous damping coefficient (normal direction) equal to $\\frac{c_{n}}{C_{n,crit}}$."))
((Real,betas,0.0,,"Fraction of the viscous damping coefficient (shear direction) equal to $\\frac{c_{s}}{C_{s,crit}}$."))
+ // Contact damping coefficient for non-linear elastic contact law (of Hertz-Mindlin type)
+ ((Real,alpha,0.0,,"Constant coefficient to define contact viscous damping for non-linear elastic force-displacement relationship."))
// temporary
((Vector3r,prevU,Vector3r::Zero(),,"Previous local displacement; only used with :yref:`Law2_L3Geom_FrictPhys_HertzMindlin`."))
((Vector2r,Fs,Vector2r::Zero(),,"Shear force in local axes (computed incrementally)"))
-
- //((Real,shearEnergy,0.0,,"Shear elastic potential energy"))
- //((Real,frictionDissipation,0.0,,"Energy dissipation due to sliding"))
- //((Real,normDampDissip,0.0,,"Energy dissipation due to sliding"))
- //((Real,shearDampDissip,0.0,,"Energy dissipation due to sliding"))
,
createIndex());
- REGISTER_CLASS_INDEX(MindlinPhys,FrictPhysTransitory);
+ REGISTER_CLASS_INDEX(MindlinPhys,FrictPhys);
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(MindlinPhys);
@@ -65,8 +69,10 @@
public :
virtual void go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction);
FUNCTOR2D(FrictMat,FrictMat);
- YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictMat_MindlinPhys,IPhysFunctor,"Calculate some physical parameters needed to obtain the normal and shear stiffnesses according to the Hertz-Mindlin's formulation (as implemented in PFC).\n\nViscous parameters can be specified either using coefficients of restitution ($e_n$, $e_s$) or viscous damping coefficient ($\\beta_n$, $\\beta_s$). The following rules apply:\n\n#. If the $\\beta_n$ ($\\beta_s$) coefficient is given, it is assigned to :yref:`MindlinPhys.betan` (:yref:`MindlinPhys.betas`) directly.\n#. If $e_n$ is given, :yref:`MindlinPhys.betan` is computed using $\\beta_n=-(\\log e_n)/\\sqrt{\\pi^2+(\\log e_n)^2}$. The same applies to $e_s$, :yref:`MindlinPhys.betas`.\n#. It is an error (exception) to specify both $e_n$ and $\\beta_n$ ($e_s$ and $\\beta_s$).\n#. If neither $e_n$ nor $\\beta_n$ is given, zero value for :yref:`MindlinPhys.betan` is used; there will be no viscous effects.\n#. If neither $e_s$ nor $\\beta_s$ is given, the value of :yref:`MindlinPhys.betan` is used for :yref:`MindlinPhys.betas` as well.\n\nThe $e_n$, $\\beta_n$, $e_s$, $\\beta_s$ are :yref:`MatchMaker` objects; they can be constructed from float values to always return constant value.\n\nSee :ysrc:`scripts/shots.py` for an example of specifying $e_n$ based on combination of parameters.",
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictMat_MindlinPhys,IPhysFunctor,"Calculate some physical parameters needed to obtain the normal and shear stiffnesses according to the Hertz-Mindlin's formulation (as implemented in PFC).\n\nViscous parameters can be specified either using coefficients of restitution ($e_n$, $e_s$) or viscous damping coefficient ($\\beta_n$, $\\beta_s$). The following rules apply:\n#. If the $\\beta_n$ ($\\beta_s$) coefficient is given, it is assigned to :yref:`MindlinPhys.betan` (:yref:`MindlinPhys.betas`) directly.\n#. If $e_n$ is given, :yref:`MindlinPhys.betan` is computed using $\\beta_n=-(\\log e_n)/\\sqrt{\\pi^2+(\\log e_n)^2}$. The same applies to $e_s$, :yref:`MindlinPhys.betas`.\n#. It is an error (exception) to specify both $e_n$ and $\\beta_n$ ($e_s$ and $\\beta_s$).\n#. If neither $e_n$ nor $\\beta_n$ is given, zero value for :yref:`MindlinPhys.betan` is used; there will be no viscous effects.\n#.If neither $e_s$ nor $\\beta_s$ is given, the value of :yref:`MindlinPhys.betan` is used for :yref:`MindlinPhys.betas` as well.\n\nThe $e_n$, $\\beta_n$, $e_s$, $\\beta_s$ are :yref:`MatchMaker` objects; they can be constructed from float values to always return constant value.\n\nSee :ysrc:`scripts/test/shots.py` for an example of specifying $e_n$ based on combination of parameters.",
((Real,gamma,0.0,,"Surface energy parameter [J/m^2] per each unit contact surface, to derive DMT formulation from HM"))
+ ((Real,eta,0.0,,"Coefficient to determinate the plastic bending moment"))
+ ((Real,krot,0.0,,"Rotational stiffness for moment contact law"))
((shared_ptr<MatchMaker>,en,,,"Normal coefficient of restitution $e_n$."))
((shared_ptr<MatchMaker>,es,,,"Shear coefficient of restitution $e_s$."))
((shared_ptr<MatchMaker>,betan,,,"Normal viscous damping coefficient $\\beta_n$."))
@@ -76,6 +82,7 @@
};
REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_MindlinPhys);
+
class Law2_ScGeom_MindlinPhys_MindlinDeresiewitz: public LawFunctor{
public:
virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
@@ -97,46 +104,46 @@
};
REGISTER_SERIALIZABLE(Law2_ScGeom_MindlinPhys_HertzWithLinearShear);
+
/******************** Law2_ScGeom_MindlinPhys_Mindlin *********/
class Law2_ScGeom_MindlinPhys_Mindlin: public LawFunctor{
// just for the deprecation warning, can be removed later
- Real _beta_parameters_of_Ip2_FrictMat_FrictMat_MindlinPhys;
+ Real _beta_parameters_of_Ip2_FrictMat_FrictMat_MindlinPhys;
bool _nothing;
public:
virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
Real normElastEnergy();
+ Real adhesionEnergy();
+
Real getfrictionDissipation();
Real getshearEnergy();
Real getnormDampDissip();
Real getshearDampDissip();
Real contactsAdhesive();
-
- FUNCTOR2D(ScGeom,MindlinPhys);
+ FUNCTOR2D(ScGeom6D,MindlinPhys);
YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(Law2_ScGeom_MindlinPhys_Mindlin,LawFunctor,"Constitutive law for the Hertz-Mindlin formulation. It includes non linear elasticity in the normal direction as predicted by Hertz for two non-conforming elastic contact bodies. In the shear direction, instead, it reseambles the simplified case without slip discussed in Mindlin's paper, where a linear relationship between shear force and tangential displacement is provided. Finally, the Mohr-Coulomb criterion is employed to established the maximum friction force which can be developed at the contact. Moreover, it is also possible to include the effect of linear viscous damping through the definition of the parameters $\\beta_{n}$ and $\\beta_{s}$.",
((bool,preventGranularRatcheting,true,,"bool to avoid granular ratcheting"))
((bool,includeAdhesion,false,,"bool to include the adhesion force following the DMT formulation. If true, also the normal elastic energy takes into account the adhesion effect."))
((bool,calcEnergy,false,,"bool to calculate energy terms (shear potential energy, dissipation of energy due to friction and dissipation of energy due to normal and tangential damping)"))
- // FIXME: all the energy attributes should be openMPAccumulator
- ((OpenMPAccumulator<Real>,frictionDissipation,,Attr::noSave,"Energy dissipation due to sliding"))
+ ((bool,includeMoment,false,,"bool to consider the rolling resistance"))
+ ((bool,LinDamp,true,,"bool to activate linear viscous damping (if false, en and es have to be defined in place of betan and betas)"))
+
+ ((OpenMPAccumulator<Real>,frictionDissipation,,Attr::noSave,"Energy dissipation due to sliding"))
((OpenMPAccumulator<Real>,shearEnergy,,Attr::noSave,"Shear elastic potential energy"))
- ((OpenMPAccumulator<Real>,normDampDissip,,Attr::noSave,"Energy dissipation due to sliding"))
- ((OpenMPAccumulator<Real>,shearDampDissip,,Attr::noSave,"Energy dissipation due to sliding"))
+ ((OpenMPAccumulator<Real>,normDampDissip,,Attr::noSave,"Energy dissipated by normal damping"))
+ ((OpenMPAccumulator<Real>,shearDampDissip,,Attr::noSave,"Energy dissipated by tangential damping"))
+
, /*deprec*/
((betan,_beta_parameters_of_Ip2_FrictMat_FrictMat_MindlinPhys,"!Moved to MindlinPhys, where the value is assigned by the appropriate Ip2 functor."))
((betas,_beta_parameters_of_Ip2_FrictMat_FrictMat_MindlinPhys,"!Moved to MindlinPhys, where the value is assigned by the appropriate Ip2 functor."))
- ((useDamping,_nothing,"Damping is now turned on automatically if either of MindlinPhys.betan or MindlinPhys.betas are non-zero."))
+ ((useDamping,_nothing,"Damping is now turned on automatically if either of MindlinPhys.betan or MindlinPhys.betas or MindlinPhys.alpha are non-zero."))
, /* init */
, /* ctor */
, /* py */
.def("contactsAdhesive",&Law2_ScGeom_MindlinPhys_Mindlin::contactsAdhesive,"Compute total number of adhesive contacts.")
-
- .def("normElastEnergy",&Law2_ScGeom_MindlinPhys_Mindlin::normElastEnergy,"Compute normal elastic potential energy. It handle the DMT formulation if :yref:`Law2_ScGeom_MindlinPhys_Mindlin::includeAdhesion` is set to true.")
- //.def("frictionDissipation",&Law2_ScGeom_MindlinPhys_Mindlin::getfrictionDissipation,"Energy dissipated by frictional behavior.")
- //.def("shearEnergy",&Law2_ScGeom_MindlinPhys_Mindlin::getshearEnergy,"Shear elastic potential energy.")
- //.def("normDampDissip",&Law2_ScGeom_MindlinPhys_Mindlin::getnormDampDissip,"Energy dissipated by normal damping.")
- //.def("shearDampDissip",&Law2_ScGeom_MindlinPhys_Mindlin::getshearDampDissip,"Energy dissipated by tangential damping.")
+ .def("normElastEnergy",&Law2_ScGeom_MindlinPhys_Mindlin::normElastEnergy,"Compute normal elastic potential energy. It handles the DMT formulation if :yref:`Law2_ScGeom_MindlinPhys_Mindlin::includeAdhesion` is set to true.")
);
DECLARE_LOGGER;
};
