yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06708
[Branch ~yade-dev/yade/trunk] Rev 2653: - Add incremental formulation for the moment rotation (bending only) to HM law (it is very simila...
------------------------------------------------------------
revno: 2653
committer: Chiara Modenese <c.modenese@xxxxxxxxx>
branch nick: yade
timestamp: Thu 2011-01-13 16:46:20 +0000
message:
- Add incremental formulation for the moment rotation (bending only) to HM law (it is very similar to the code for the shear part as written in ScGeom - it has been tested but feedbacks are welcome).
PS Let me know if you are interested to have it also in "Law2_ScGeom6D_CohFrictPhys_CohesionMoment" where actually the total formulation is already present.
@Bruno: maybe in ScGeom we could add this incremental formulation too?
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-11 19:01:11 +0000
+++ pkg/dem/HertzMindlin.cpp 2011-01-13 16:46:20 +0000
@@ -70,6 +70,7 @@
/* pass values calculated from above to MindlinPhys */
mindlinPhys->tangensOfFrictionAngle = std::tan(frictionAngle);
+ //mindlinPhys->prevNormal = scg->normal; // used to compute relative rotation
mindlinPhys->kno = Kno; // this is just a coeff
mindlinPhys->kso = Kso; // this is just a coeff
mindlinPhys->adhesionForce = Adhesion;
@@ -112,7 +113,7 @@
Real normEnergy=0;
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
if(!I->isReal()) continue;
- ScGeom6D* scg = dynamic_cast<ScGeom6D*>(I->geom.get());
+ ScGeom* scg = dynamic_cast<ScGeom*>(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);}
@@ -128,7 +129,7 @@
Real adhesionEnergy=0;
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
if(!I->isReal()) continue;
- ScGeom6D* scg = dynamic_cast<ScGeom6D*>(I->geom.get());
+ ScGeom* scg = dynamic_cast<ScGeom*>(I->geom.get());
MindlinPhys* phys = dynamic_cast<MindlinPhys*>(I->phys.get());
if (phys && includeAdhesion) {
Real R=scg->radius1*scg->radius2/(scg->radius1+scg->radius2);
@@ -230,7 +231,7 @@
State* de1 = Body::byId(id1,scene)->state.get();
State* de2 = Body::byId(id2,scene)->state.get();
- ScGeom6D* scg = static_cast<ScGeom6D*>(ig.get());
+ ScGeom* scg = static_cast<ScGeom*>(ig.get());
MindlinPhys* phys = static_cast<MindlinPhys*>(ip.get());
const shared_ptr<Body>& b1=Body::byId(id1,scene);
@@ -412,19 +413,63 @@
/********************************************/
/* MOMENT CONTACT LAW (only bending moment) */
/********************************************/
-
- if (includeMoment){
- phys->moment_bending = scg->getBending()*phys->kr;
+ if (includeMoment){
+ // new code to compute relative particle rotation (similar to the way the shear is computed)
+ Vector3r relAngVel = (b2->state->angVel-b1->state->angVel);
+ relAngVel = relAngVel - scg->normal.dot(relAngVel)*scg->normal; // keep only the bending part
+ Vector3r relRot = relAngVel*dt; // relative rotation due to rolling behaviour
+ // incremental formulation for the bending moment (as for the shear part)
+ Vector3r& momentBend = phys->momentBend;
+ momentBend = scg->rotate(momentBend); // rotate moment vector (updated)
+ momentBend = momentBend-phys->kr*relRot; // add incremental rolling to the rolling vector FIXME: is the sign correct?
+
+#if 0
+ // code to compute the relative particle rotation
+ if (includeMoment){
+ Real rMean = (scg->radius1+scg->radius2)/2.;
+ // sliding motion
+ Vector3r duS1 = scg->radius1*(phys->prevNormal-scg->normal);
+ Vector3r duS2 = scg->radius2*(scg->normal-phys->prevNormal);
+ // rolling motion
+ Vector3r duR1 = scg->radius1*dt*b1->state->angVel.cross(scg->normal);
+ Vector3r duR2 = -scg->radius2*dt*b2->state->angVel.cross(scg->normal);
+ // relative position of the old contact point with respect to the new one
+ Vector3r relPosC1 = duS1+duR1;
+ Vector3r relPosC2 = duS2+duR2;
+
+ Vector3r duR = (relPosC1+relPosC2)/2.; // incremental displacement vector (same radius is temporarily assumed)
+
+ // check wheter rolling will be present, if not do nothing
+ Vector3r x=scg->normal.cross(duR);
+ Vector3r normdThetaR(Vector3r::Zero()); // initialize
+ if(x.squaredNorm()==0) { /* no rolling */ }
+ else {
+ Vector3r normdThetaR = x/x.norm(); // moment unit vector
+ phys->dThetaR = duR.norm()/rMean*normdThetaR;} // incremental rolling
+
+ // incremental formulation for the bending moment (as for the shear part)
+ Vector3r& momentBend = phys->momentBend;
+ momentBend = scg->rotate(momentBend); // rotate moment vector
+ momentBend = momentBend+phys->kr*phys->dThetaR; // add incremental rolling to the rolling vector FIXME: is the sign correct?
+#endif
+
+ // check plasticity condition
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;
+ Real scalarMoment = phys->momentBend.norm();
+ if (MomentMax>0){
+ if(scalarMoment > MomentMax)
+ {
+ Real ratio = MomentMax/scalarMoment; // to fix the moment to its yielding value
+ phys->momentBend *= ratio;
+ }
}
- scene->forces.addTorque(id1,-phys->moment_bending);
- scene->forces.addTorque(id2, phys->moment_bending);
+ // apply moments
+ scene->forces.addTorque(id1,-phys->momentBend);
+ scene->forces.addTorque(id2,phys->momentBend);
}
+
+ // update variables
+ //phys->prevNormal = scg->normal;
}
=== modified file 'pkg/dem/HertzMindlin.hpp'
--- pkg/dem/HertzMindlin.hpp 2011-01-11 19:01:11 +0000
+++ pkg/dem/HertzMindlin.hpp 2011-01-13 16:46:20 +0000
@@ -39,7 +39,10 @@
((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"))
+
+ //((Vector3r,prevNormal,Vector3r::Zero(),,"Save previous contact normal to compute relative rotation"))
+ ((Vector3r,momentBend,Vector3r::Zero(),,"Artificial bending moment to provide rolling resistance in order to account for some degree of interlocking between particles"))
+ //((Vector3r,dThetaR,Vector3r::Zero(),,"Incremental rolling vector"))
((Real,radius,NaN,,"Contact radius (only computed with :yref:`Law2_ScGeom_MindlinPhys_Mindlin::calcEnergy`)"))
@@ -53,6 +56,7 @@
// 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)"))
@@ -64,6 +68,7 @@
REGISTER_SERIALIZABLE(MindlinPhys);
+
/******************** Ip2_FrictMat_FrictMat_MindlinPhys *******/
class Ip2_FrictMat_FrictMat_MindlinPhys: public IPhysFunctor{
public :
@@ -71,7 +76,7 @@
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#. 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,eta,0.0,,"Coefficient to determine 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$."))
@@ -122,12 +127,12 @@
Real getshearDampDissip();
Real contactsAdhesive();
- FUNCTOR2D(ScGeom6D,MindlinPhys);
+ FUNCTOR2D(ScGeom,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)"))
- ((bool,includeMoment,false,,"bool to consider the rolling resistance"))
+ ((bool,includeMoment,false,,"bool to consider rolling resistance (if :yref:`Ip2_FrictMat_FrictMat_MindlinPhys::eta` is 0.0, no plastic condition is applied.)"))
((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"))
Follow ups
