yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06790
[Branch ~yade-dev/yade/trunk] Rev 2676: Changes in "normalInelasticity" files suggested by http://www.mail-archive.com/yade-dev@xxxxxxxxx...
------------------------------------------------------------
revno: 2676
committer: jduriez <jduriez@c1solimara-l>
branch nick: yade
timestamp: Mon 2011-01-24 15:45:35 +0100
message:
Changes in "normalInelasticity" files suggested by http://www.mail-archive.com/yade-dev@xxxxxxxxxxxxxxxxxxx/msg06011.html
Add of some words in the doc of normal in pkg/dem/DemXDofGeom.hpp
modified:
doc/yade-articles.bib
pkg/dem/DemXDofGeom.hpp
pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp
pkg/dem/KinemCNSEngine.cpp
pkg/dem/NormalInelasticityLaw.cpp
pkg/dem/NormalInelasticityPhys.hpp
scripts/normalInelasticity-test.py
--
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 'doc/yade-articles.bib'
--- doc/yade-articles.bib 2010-12-15 15:22:47 +0000
+++ doc/yade-articles.bib 2011-01-24 14:45:35 +0000
@@ -277,6 +277,7 @@
title = "A discrete modeling-based constitutive relation for infilled rock joints",
year = "2010",
journal = "International Journal of Rock Mechanics & Mining Sciences",
+ doi = "10.1016/j.ijrmms.2010.09.008"
note = "in press"
}
=== modified file 'pkg/dem/DemXDofGeom.hpp'
--- pkg/dem/DemXDofGeom.hpp 2010-12-13 12:11:43 +0000
+++ pkg/dem/DemXDofGeom.hpp 2011-01-24 14:45:35 +0000
@@ -8,7 +8,7 @@
class GenericSpheresContact: public IGeom{
YADE_CLASS_BASE_DOC_ATTRS_CTOR(GenericSpheresContact,IGeom,
"Class uniting :yref:`ScGeom` and :yref:`Dem3DofGeom`, for the purposes of :yref:`GlobalStiffnessTimeStepper`. (It might be removed inthe future). Do not use this class directly.",
- ((Vector3r,normal,,,"Unit vector oriented along the interaction. |yupdate|"))
+ ((Vector3r,normal,,,"Unit vector oriented along the interaction, from particle #1, towards particle #2. |yupdate|"))
((Vector3r,contactPoint,,,"some reference point for the interaction (usually in the middle). |ycomp|"))
((Real,refR1,,,"Reference radius of particle #1. |ycomp|"))
((Real,refR2,,,"Reference radius of particle #2. |ycomp|")),
=== modified file 'pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp'
--- pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp 2010-12-31 14:35:21 +0000
+++ pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp 2011-01-24 14:45:35 +0000
@@ -44,30 +44,26 @@
Real Eb = sdec2->young;
Real Va = sdec1->poisson;
Real Vb = sdec2->poisson;
- Real Da = geom->radius1; // FIXME - multiply by factor of sphere interaction distance (so sphere interacts at bigger range that its geometrical size)
- Real Db = geom->radius2; // FIXME - as above
+ Real Ra = geom->radius1; // FIXME - multiply by factor of sphere interaction distance (so sphere interacts at bigger range that its geometrical size)
+ Real Rb = geom->radius2; // FIXME - as above
Real fa = sdec1->frictionAngle;
Real fb = sdec2->frictionAngle;
- Real Kn = 2.0*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
-
- Real Ks = 2.0*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);//harmonic average of two stiffnesses with ks=V*kn for each sphere
-
- // Jean-Patrick Plassiard, Noura Belhaine, Frederic
-// Victor Donze, "Calibration procedure for spherical
- // discrete elements using a local moemnt law".
- Real Kr = betaR*std::pow((Da+Db)/2.0,2)*Ks;
-
- contactPhysics->frictionAngle = std::min(fa,fb); // FIXME - this is actually a waste of memory space
- contactPhysics->tangensOfFrictionAngle = std::tan(contactPhysics->frictionAngle);
- contactPhysics->forMaxMoment = 1.0*(Da+Db)/2.0; // 1.0 corresponding to ethaR which I don't know exactly where to define as a parameter...
+ Real Kn = 2.0*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);//harmonic average of two stiffnesses
+
+ Real Ks = 2.0*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Va);//harmonic average of two stiffnesses with ks=V*kn for each sphere
+
+ // Jean-Patrick Plassiard, Noura Belheine, Frederic Victor Donze, "A Spherical Discrete Element Model: calibration procedure and incremental response", DOI: 10.1007/s10035-009-0130-x
+
+ Real Kr = betaR*std::pow((Ra+Rb)/2.0,2)*Ks;
+
+ contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb));
+ contactPhysics->forMaxMoment = 1.0*(Ra+Rb)/2.0; // 1.0 corresponding to ethaR which I don't know exactly where to define as a parameter...
// Lot of suppress here around (>) r2276. Normally not bad but ? See Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp to re-find the initial §...
-
- contactPhysics->prevNormal = geom->normal;
-
+
contactPhysics->knLower = Kn;
- contactPhysics->kn = Kn;
+ contactPhysics->kn = Kn;
contactPhysics->ks = Ks;
contactPhysics->kr = Kr;
}
=== modified file 'pkg/dem/KinemCNSEngine.cpp'
--- pkg/dem/KinemCNSEngine.cpp 2010-11-07 11:46:20 +0000
+++ pkg/dem/KinemCNSEngine.cpp 2011-01-24 14:45:35 +0000
@@ -44,121 +44,6 @@
}
-/*void KinemCNSEngine::letMove()
-{
- shared_ptr<BodyContainer> bodies = scene->bodies;
-
- if(LOG) cout << "It : " << scene->iter << endl;
-
-
- const Real& dt = scene->dt;
- Real dx = shearSpeed * dt;
-
-
- Real Ysup = topbox->state->pos.y();
- Real Ylat = leftbox->state->pos.y();
-
-// Changes in vertical and horizontal position :
-
- topbox->state->pos += Vector3r(dx,deltaH,0);
-
- leftbox->state->pos += Vector3r(dx/2.0,deltaH/2.0,0);
- rightbox->state->pos += Vector3r(dx/2.0,deltaH/2.0,0);
-
- Real Ysup_mod = topbox->state->pos.y();
- Real Ylat_mod = leftbox->state->pos.y();
-
-// with the corresponding velocities :
- topbox->state->vel = Vector3r(shearSpeed,deltaH/dt,0);
- leftbox->state->vel = Vector3r(shearSpeed/2.0,deltaH/(2.0*dt),0);
- rightbox->state->vel = Vector3r(shearSpeed/2.0,deltaH/(2.0*dt),0);
-
-// Then computation of the angle of the rotation to be done :
- if (alpha == Mathr::PI/2.0) // Case of the very beginning
- {
- dalpha = - atan( dx / (Ysup_mod -Ylat_mod) );
- }
- else
- {
- Real A = (Ysup_mod - Ylat_mod) * 2.0*tan(alpha) / (2.0*(Ysup - Ylat) + dx*tan(alpha) );
- dalpha = atan( (A - tan(alpha))/(1.0 + A * tan(alpha)));
- }
-
- Quaternionr qcorr(AngleAxisr(dalpha,Vector3r::UnitZ()));
- if(LOG)
- cout << "Quaternion associe a la rotation incrementale : " << qcorr.w() << " " << qcorr.x() << " " << qcorr.y() << " " << qcorr.z() << endl;
-// On applique la rotation en changeant l'orientation des plaques, leurs vang et en affectant donc alpha
- leftbox->state->ori = qcorr*leftbox->state->ori;
- leftbox->state->angVel = Vector3r(0,0,1)*dalpha/dt;
-
- rightbox->state->ori = qcorr*rightbox->state->ori;
- rightbox->state->angVel = Vector3r(0,0,1)*dalpha/dt;
-
- gamma+=dx;
-}*/
-
-
-// void KinemCNSEngine::computeDu()
-// {
-//
-// scene->forces.sync(); Vector3r F_sup=scene->forces.getForce(id_topbox);
-//
-// if(firstRun)
-// {
-// if ( !myLdc ) // FIXME : reenable this feature since Law2 is a functor
-// {
-// // vector<shared_ptr<Engine> >::iterator itFirst = scene->engines.begin();
-// // vector<shared_ptr<Engine> >::iterator itLast = scene->engines.end();
-// // for ( ;itFirst!=itLast; ++itFirst )
-// // {
-// // if ( ( *itFirst )->getClassName() == "Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity" )
-// // {
-// // myLdc = YADE_PTR_CAST<Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity> ( *itFirst );
-// // coeff_dech = myLdc ->coeff_dech;
-// // if(LOG) cout << "Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity engine found, with coeff_dech = " << coeff_dech << endl;
-// // }
-// // }
-// }
-//
-// alpha=Mathr::PI/2.0;;
-// Y0 = topbox->state->pos.y();
-// cout << "Y0 initialise à : " << Y0 << endl;
-// F_0 = F_sup.y();
-// firstRun=false;
-// }
-//
-//
-// computeStiffness();
-// Real Hcurrent = topbox->state->pos.y();
-// Real Fdesired = F_0 + KnC * 1.0e9 * Scontact * (Hcurrent-Y0); // The value of the force desired
-//
-// // Prise en compte de la difference de rigidite entre charge et decharge dans le cadre de Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity : => INUTILE maintenant ?
-// // if( F_sup.y() > Fdesired ) // cas ou l'on va monter la plaq <=> (normalemt) a une decharge
-// // stiffness *= coeff_dech;
-//
-// if( (KnC==0) && (stiffness==0) )
-// {
-// deltaH=0;
-// cerr << "KnC et stiffness(sample) = 0 => DNC en fait et non CNC..." << endl;
-// }
-// else
-// {
-// deltaH = ( F_sup.y() - ( Fdesired ))/(stiffness+KnC* 1.0e9 * Scontact);
-// }
-//
-// if(LOG) cout << "Alors q je veux KnC = " << KnC << " depuis F_0 = " << F_0 << " et Y0 = " << Y0 << endl;
-// if(LOG) cout << "deltaH a permettre normalement :" << deltaH << endl;
-//
-// deltaH = (1-wallDamping)*deltaH;
-// if(LOG) cout << "deltaH apres amortissement :" << deltaH << endl;
-//
-// if(abs(deltaH) > max_vel*scene->dt)
-// {
-// deltaH=deltaH/abs(deltaH)*max_vel*scene->dt;
-// if(LOG) cout << "Correction appliquee pour ne pas depasser vmax(comp)" << endl;
-// }
-//
-// }
=== modified file 'pkg/dem/NormalInelasticityLaw.cpp'
--- pkg/dem/NormalInelasticityLaw.cpp 2010-12-31 14:35:21 +0000
+++ pkg/dem/NormalInelasticityLaw.cpp 2011-01-24 14:45:35 +0000
@@ -18,17 +18,12 @@
void Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity::go(shared_ptr<IGeom>& iG, shared_ptr<IPhys>& iP, Interaction* contact)
{
-
- const Real& dt = scene->dt;
-
int id1 = contact->getId1();
int id2 = contact->getId2();
// cout << "contact entre " << id1 << " et " << id2;
- State* de1 = Body::byId(id1,scene)->state.get();
- State* de2 = Body::byId(id2,scene)->state.get();
NormalInelasticMat* Mat1 = static_cast<NormalInelasticMat*>(Body::byId(id1,scene)->material.get());
- ScGeom6D* currentContactGeometry = YADE_CAST<ScGeom6D*>(iG.get());
+ ScGeom6D* geom = YADE_CAST<ScGeom6D*>(iG.get());
NormalInelasticityPhys* currentContactPhysics = YADE_CAST<NormalInelasticityPhys*> (iP.get());
Vector3r& shearForce = currentContactPhysics->shearForce;
@@ -44,9 +39,9 @@
// ******** Computation of normal Force : depends of the history ******* //
- Real Fn; // la valeur de Fn qui va etre calculee selon différentes manieres puis affectee
+ Real Fn; // Fn's value, computed by different means
// cout << " Dans Law2 valeur de kn : " << currentContactPhysics->kn << endl;
- Real un = currentContactGeometry->penetrationDepth; // >0 for real penetration
+ Real un = geom->penetrationDepth; // >0 for real penetration
// cout << "un = " << un << " alors que unMax = "<< currentContactPhysics->unMax << " et previousun = " << currentContactPhysics->previousun << " et previousFn =" << currentContactPhysics->previousFn << endl;
@@ -73,7 +68,7 @@
// cout << "j'ai corrige pour ne pas etre negatif" << endl;
}
}
- currentContactPhysics->normalForce = Fn*currentContactGeometry->normal;
+ currentContactPhysics->normalForce = Fn*geom->normal;
// cout << "Fn appliquee " << Fn << endl << endl;
// actualisation :
@@ -94,48 +89,10 @@
}
else
{
- Vector3r axis;
- Real angle;
-
-// Here is the code with approximated rotations
-
- axis = currentContactPhysics->prevNormal.cross(currentContactGeometry->normal);
- shearForce -= shearForce.cross(axis);
- angle = dt*0.5*currentContactGeometry->normal.dot(de1->angVel+de2->angVel);
- axis = angle*currentContactGeometry->normal;
- shearForce -= shearForce.cross(axis);
-
-// Here is the code with exact rotationns
-
- // Quaternionr q;
- //
- // axis = currentContactPhysics->prevNormal.cross(currentContactGeometry->normal);
- // angle = acos(currentContactGeometry->normal.dot(currentContactPhysics->prevNormal));
- // q.fromAngleAxis(angle,axis);
- //
- // currentContactPhysics->shearForce = currentContactPhysics->shearForce*q;
- //
- // angle = dt*0.5*currentContactGeometry->normal.dot(de1->angVel+de2->angVel);
- // axis = currentContactGeometry->normal;
- // q.fromAngleAxis(angle,axis);
- // currentContactPhysics->shearForce = q*currentContactPhysics->shearForce;
-
-
- Vector3r x = currentContactGeometry->contactPoint;
- Vector3r c1x = (x - de1->se3.position);
- Vector3r c2x = (x - de2->se3.position);
- /// 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, or discussions on yade mailing lists)
- Vector3r _c1x_ = currentContactGeometry->radius1*currentContactGeometry->normal;
- Vector3r _c2x_ = -currentContactGeometry->radius2*currentContactGeometry->normal;
- Vector3r relativeVelocity = (de2->vel+de2->angVel.cross(_c2x_)) - (de1->vel+de1->angVel.cross(_c1x_));
- Vector3r shearVelocity = relativeVelocity-currentContactGeometry->normal.dot(relativeVelocity)*currentContactGeometry->normal;
- Vector3r shearDisplacement = shearVelocity*dt;
-
-
- shearForce -= currentContactPhysics->ks*shearDisplacement;
+ // Correction of previous shear force to take into account the change in normal:
+ shearForce = geom->rotate(currentContactPhysics->shearForce);
+ // Update of shear force corresponding to shear displacement increment:
+ shearForce -= currentContactPhysics->ks * geom->shearIncrement();
Real maxFs = 0;
Real Fs = currentContactPhysics->shearForce.norm();
@@ -156,8 +113,8 @@
Vector3r f = currentContactPhysics->normalForce + shearForce;
scene->forces.addForce (id1,-f);
scene->forces.addForce (id2, f);
- scene->forces.addTorque(id1,-c1x.cross(f));
- scene->forces.addTorque(id2, c2x.cross(f));
+ scene->forces.addTorque(id1,-(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(f));
+ scene->forces.addTorque(id2, -(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(f));
// *** End of computation of tangential force *** //
@@ -165,8 +122,8 @@
if(momentRotationLaw)
{
- currentContactPhysics->moment_twist = (currentContactGeometry->getTwist()*currentContactPhysics->kr)*currentContactGeometry->normal ;
- currentContactPhysics->moment_bending = currentContactGeometry->getBending() * currentContactPhysics->kr;
+ currentContactPhysics->moment_twist = (geom->getTwist()*currentContactPhysics->kr)*geom->normal ;
+ currentContactPhysics->moment_bending = geom->getBending() * currentContactPhysics->kr;
moment = currentContactPhysics->moment_twist + currentContactPhysics->moment_bending;
@@ -184,8 +141,6 @@
scene->forces.addTorque(id2, moment);
}
// ******** Moment law END ******* //
-
- currentContactPhysics->prevNormal = currentContactGeometry->normal;
}
}
=== modified file 'pkg/dem/NormalInelasticityPhys.hpp'
--- pkg/dem/NormalInelasticityPhys.hpp 2010-12-31 14:35:21 +0000
+++ pkg/dem/NormalInelasticityPhys.hpp 2011-01-24 14:45:35 +0000
@@ -11,12 +11,12 @@
#include<yade/pkg/dem/FrictPhys.hpp>
-class NormalInelasticityPhys : public FrictPhysTransitory
+class NormalInelasticityPhys : public FrictPhys
{
public :
virtual ~NormalInelasticityPhys();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(NormalInelasticityPhys,FrictPhysTransitory,
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(NormalInelasticityPhys,FrictPhys,
"Physics (of interaction) for using :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity` : with inelastic unloadings",
((Real,unMax,0.0,,"the maximum value of penetration depth of the history of this interaction"))
((Real,previousun,0.0,,"the value of this un at the last time step"))
@@ -30,7 +30,7 @@
,
createIndex();
);
- REGISTER_CLASS_INDEX(NormalInelasticityPhys,FrictPhysTransitory);
+ REGISTER_CLASS_INDEX(NormalInelasticityPhys,FrictPhys);
};
REGISTER_SERIALIZABLE(NormalInelasticityPhys);
=== modified file 'scripts/normalInelasticity-test.py'
--- scripts/normalInelasticity-test.py 2010-11-08 15:26:46 +0000
+++ scripts/normalInelasticity-test.py 2011-01-24 14:45:35 +0000
@@ -7,7 +7,7 @@
#Different graphs illustrate the effects of the different loadings. The run is paused at each plot window (so that there is time to observe them). Push on "Return", while being in the Yade terminal, to go ahead.
-#All comments concern r2528
+#No crash warranty with r2676
from yade import plot
@@ -16,8 +16,8 @@
O.materials.append(NormalInelasticMat(density=2600,young=4.0e9,poisson=.04,frictionAngle=.6,coeff_dech=3.0,label='Materiau1'))
#Def of the bodies of the simulations : 2 spheres, with names which will be useful after
-O.bodies.append(utils.sphere([0,0,0], 1, dynamic=False, wire=False, color=None, highlight=False)) #'Materiau1', as the latest material defined, will be used
-O.bodies.append(utils.sphere([0,2,0], 1, dynamic=False, wire=False, color=None, highlight=False))
+O.bodies.append(utils.sphere([0,0,0], 1, fixed=True, wire=False, color=None, highlight=False)) #'Materiau1', as the latest material defined, will be used
+O.bodies.append(utils.sphere([0,2,0], 1, fixed=True, wire=False, color=None, highlight=False))
lowerSphere=O.bodies[0]
upperSphere=O.bodies[1]
@@ -26,7 +26,7 @@
#Def of the engines taking part to the simulation loop
O.engines=[
ForceResetter(),
- InsertionSortCollider([Bo1_Sphere_Aabb()]),
+ InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=0),# use of verletDist>0 without NewtonIntegrator let crash ??
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D()],
[Ip2_2xNormalInelasticMat_NormalInelasticityPhys(betaR=0.24)],
@@ -37,7 +37,7 @@
#Def of the python command letMove() :
-# Will move "by hand" the upperSphere towards or away from the lower one. Modifying by hand only the speed of bodies is indeed not sufficient, see NewtonsIntegrator, and https://bugs.launchpad.net/yade/+bug/398089. Alternative way is presented below
+# Will move "by hand" the upperSphere towards or away from the lower one.
def letMove():#Load for the first 10 iterations, unload for the 7 following iterations, then reload
vImposed=[0,-1,0]
if O.iter < 25 and O.iter>14:
@@ -107,7 +107,7 @@
# - un decreases indeed during this shear, but maybe a zoom on the curves is needed to see it.
# - We can observe that the force state of the sample decreases a line with a slope equal to tan(~34.5°)=tan(~0.602 rad). Why not strict equality ? Because of the measure of the slope or because something else ? To see...
# - during this phase O.forces.t(0).norm() / O.forces.f(0)[0] seems to increase between 0.502 and 0.507 (according to r=[T[i]/F[i] for i in range(50,T.__len__()-20) ])
-# Note to explain this that Fx = O.forces.f(0)[0] is more and more different from Ft, from which we can expect Ft = Torque /Radius
+# Please note, to explain this, that Fx = O.forces.f(0)[0] is more and more different from Ft, from which we can expect Ft = Torque /Radius
## ------ Test of the law for the moment, using blockedDOF_s ------ #
#O.loadTmp('EndComp')
