yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #04013
[Branch ~yade-dev/yade/trunk] Rev 2146: - Add of a script NormalInelasticityTest to test ... the NormalInelasticityLaw (Law2_ScGeom_Norma...
------------------------------------------------------------
revno: 2146
committer: jduriez <jduriez@c1solimara-l>
branch nick: trunk
timestamp: Thu 2010-04-15 18:33:51 +0200
message:
- Add of a script NormalInelasticityTest to test ... the NormalInelasticityLaw (Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity)
Moreover my python problems I had difficulties to obtain "perfect results". I think it is linked to the way un is computed in Ig2_Sphere_Sphere_ScGeom (see for
example this "shift2").
Why not making penetrationDepth directly accessible through python ??
- As discussed with Chiara directly, typos in mindlin_test
added:
scripts/NormalInelasticityTest.py
modified:
pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp
pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.hpp
scripts/mindlin_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 'pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp 2010-04-08 09:14:11 +0000
+++ pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp 2010-04-15 16:33:51 +0000
@@ -21,6 +21,7 @@
void Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity::action()// a remplacer par :
// void Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity::go(shared_ptr<InteractionGeometry>& iG, shared_ptr<InteractionPhysics>& iP, Interaction* contact, Scene* scene)
{
+// cout << "\n Nvlle it :"<< endl;
shared_ptr<BodyContainer>& bodies = scene->bodies;
Real dt = Omega::instance().getTimeStep();
@@ -38,7 +39,7 @@
const shared_ptr<Interaction>& contact = *ii;// supprimable
int id1 = contact->getId1();
int id2 = contact->getId2();
- cout << "contact entre " << id1 << " et " << id2 << " reel ? " << contact->isReal() << endl;
+// cout << "contact entre " << id1 << " et " << id2 << " reel ? " << contact->isReal() << endl;
if ( !( (*bodies)[id1]->getGroupMask() & (*bodies)[id2]->getGroupMask() & sdecGroupMask) )
continue; // skip other groups,
@@ -61,6 +62,7 @@
Real previousun = currentContactPhysics->previousun;
Real previousFn = currentContactPhysics->previousFn;
Real kn = currentContactPhysics->kn;
+// cout << "Valeur de kn : " << kn << endl;
Real Fn; // la valeur de Fn qui va etre calculee selon différentes manieres puis affectee
Real un = currentContactGeometry->penetrationDepth; // compte positivement lorsq vraie interpenetration
@@ -69,17 +71,17 @@
{
Fn = kn*un;
currentContactPhysics->unMax = std::abs(un);
- cout << "je suis dans le calcul normal " << endl;
+// cout << "je suis dans le calcul normal " << endl;
}
else
{
Fn = previousFn + coeff_dech * kn * (un-previousun); // Calcul incrémental de la nvlle force
- cout << "je suis dans l'autre calcul" << endl;
+// cout << "je suis dans l'autre calcul" << endl;
if(std::abs(Fn) > std::abs(kn * un)) // verif qu'on ne depasse la courbe
Fn = kn*un;
if(Fn < 0.0 ) // verif qu'on reste positif FIXME ATTENTION ON S'EST FICHU DU NORMAL ADHESION !!!!
{Fn = 0;
- cout << "j'ai corrige pour ne pas etre negatif" << endl;
+// cout << "j'ai corrige pour ne pas etre negatif" << endl;
}
}
@@ -150,45 +152,33 @@
shearForce -= currentContactPhysics->ks*shearDisplacement;
-/////////////////////// CREEP END /////////////////////////////
- // cerr << "shearForce = " << shearForce << endl;
Real maxFs = 0;
Real Fs = currentContactPhysics->shearForce.Length();
maxFs = std::max((Real) 0,Fn*currentContactPhysics->tangensOfFrictionAngle);
if ( Fs > maxFs )
{
-
-
+
currentContactPhysics->SetBreakingState();
- // maxFs = currentContactPhysics->shearAdhesion;
- // if (!currentContactPhysics->cohesionDisablesFriction && un>0)
- // maxFs += Fn * currentContactPhysics->tangensOfFrictionAngle;
-
- maxFs = max((Real) 0, Fn * currentContactPhysics->tangensOfFrictionAngle);
-
- //cerr << "currentContactPhysics->tangensOfFrictionAngle = " << currentContactPhysics->tangensOfFrictionAngle << endl;
- // cerr << "maxFs = " << maxFs << endl;
-
- maxFs = maxFs / Fs;
- // cerr << "maxFs = " << maxFs << endl;
- if (maxFs>1)
+ maxFs = max((Real) 0, Fn * currentContactPhysics->tangensOfFrictionAngle);
+
+ maxFs = maxFs / Fs;
+ if (maxFs>1)
cerr << "maxFs>1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
shearForce *= maxFs;
if (Fn<0) currentContactPhysics->normalForce = Vector3r::ZERO;
}
-
//////// PFC3d SlipModel
- 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));
+ 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));
// Moment law ///
=== modified file 'pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.hpp 2010-04-08 09:14:11 +0000
+++ pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.hpp 2010-04-15 16:33:51 +0000
@@ -25,7 +25,7 @@
YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity,
// LawFunctor,
GlobalEngine,
- "Contact law including cohesion, moment transfer and inelastic compression behaviour\n\n This contact Law is inspired by :yref:`CohesiveFrictionalContactLaw` (inspired itselve directly from the work of Plassiard & Belheine, see the corresponding articles in (Annual Report 2006) in http://geo.hmg.inpg.fr/frederic/Discrete_Element_Group_FVD.html for example).\n\n It allows so to set moments, cohesion, tension limit and (that's the difference) inelastic unloadings in compression between bodies. All that concerned brokenBodies (this flag and the erosionactivated one) and the useless 'iter' has been suppressed.\n\n The Relationsships corresponding are Ip2_2xCohFrictMat_NormalInelasticityPhys, where the rigidities, the friction angles (with their tan()), and the orientations of the interactions are calculated. No more cohesion and tension limits are computed for all the interactions.\n\n To use it you should also use :\n- :yref:`CohesiveFrictionalMat` for bodies, with *isCohesive* = 1 (A verifier ce dernier point)\n- :yref:`Ip2_2xCohFrictMat_NormalInelasticityPhys` (=> which involves interactions of :yref:`NormalInelasticityPhys` type)",
+ "Contact law including cohesion, moment transfer and inelastic compression behaviour\n\n This contact Law is inspired by :yref:`CohesiveFrictionalContactLaw` (inspired itselve directly from the work of Plassiard & Belheine, see the corresponding articles in (Annual Report 2006) in http://geo.hmg.inpg.fr/frederic/Discrete_Element_Group_FVD.html for example).\n\n It allows so to set moments, cohesion, tension limit and (that's the difference) inelastic unloadings in compression between bodies. All that concerned brokenBodies (this flag and the erosionactivated one) and the useless 'iter' has been suppressed.\n\n The Relationsships corresponding are Ip2_2xCohFrictMat_NormalInelasticityPhys, where the rigidities, the friction angles (with their tan()), and the orientations of the interactions are calculated. No more cohesion and tension limits are computed for all the interactions.\n\n To use it you should also use :\n- :yref:`CohesiveFrictionalMat` for bodies, with *isCohesive* = 1 (A verifier ce dernier point)\n- :yref:`Ip2_2xCohFrictMat_NormalInelasticityPhys` (=> which involves interactions of :yref:`NormalInelasticityPhys` type).\n\n The effect of this law on the normal force are illustrated in scripts/NormalInelasticityTest.py",
((int,sdecGroupMask,1,"?"))
((Real,coeff_dech,1.0,"=kn(unload) / kn(load)"))
((bool,momentRotationLaw,true,"boolean, true=> computation of a torque (against relative rotation) exchanged between particles"))
=== added file 'scripts/NormalInelasticityTest.py'
--- scripts/NormalInelasticityTest.py 1970-01-01 00:00:00 +0000
+++ scripts/NormalInelasticityTest.py 2010-04-15 16:33:51 +0000
@@ -0,0 +1,61 @@
+# Script to test the constitutive law contained in NormalInelasticityLaw : consider two spheres whose penetration of the contact evolves => Monitor of the normal force
+
+from yade import plot
+
+#Def of the material which will be used
+O.materials.append(CohesiveFrictionalMat(density=2600,young=4.0e9,poisson=.04,frictionAngle=.6,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))
+
+LowerSphere=O.bodies[0]
+UpperSphere=O.bodies[1]
+
+
+#Def of the engines taking part to the simulation loop
+O.engines=[
+ ForceResetter(),
+ BoundDispatcher([Bo1_Sphere_Aabb()]),
+ InsertionSortCollider(),
+ InteractionGeometryDispatcher([Ig2_Sphere_Sphere_ScGeom()]),
+ InteractionPhysicsDispatcher([Ip2_2xCohFrictMat_NormalInelasticityPhys()]),
+ PeriodicPythonRunner(iterPeriod=1,command='letMove()'),
+ Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity(coeff_dech=3)
+]
+
+
+#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. For such purposes you could also use TranslationEngine, or you can also use JumpChangeSe3, which applies finite change in position/orientation in each step
+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:
+ vImposed=[0,1,0]
+ UpperSphere.state.vel=vImposed
+ UpperSphere.state.pos=UpperSphere.state.pos+UpperSphere.state.vel*O.dt
+
+
+#Def of the python command defData() : which will be used once the interaction will really exist
+def defData():
+ i=O.interactions[1,0]
+ VecFn=i.phys.normalForce
+ VecDist=UpperSphere.state.pos-LowerSphere.state.pos
+ plot.addData(Normfn=VecFn.Length(),FnY=VecFn[1],step=O.iter,un=LowerSphere.shape.radius+UpperSphere.shape.radius-VecDist.Length()) # the 1e-5 because of the "shift2" which appears in the computation of the penetration depth made in Ig2_Sphere_Sphere_ScGeom. It seems that this shift2 = 1e-5
+
+
+
+
+O.dt=1e-5
+
+yade.qt.View()
+O.run(2,True) #cycles "for free", so that the interaction between spheres will be defined (with his physics and so on)
+O.engines=O.engines+[PeriodicPythonRunner(iterPeriod=1,command='defData()')]
+
+#O.saveTmp('INL')
+O.run(40,True)
+
+# define of the plots to be made : un(step), and Fn(un)
+plot.plots={'step':('un',),'un':('Normfn',)}
+plot.plot()
+
+#NB : the shape of the curve Fn(un) seems to not be perfect. It is indeed not because of NormalInelasticityLaw. But of differences between the un computed here in this python script and the one which is computed in Ig2_Sphere_Sphere_ScGeom (see for example this "shift2"). Fn being linked to this last un, these slight differences explain the shape of the curves. If phys.penetrationDepth would exist in python, and thus could be directly considered, I think the curve would be perfect !
\ No newline at end of file
=== modified file 'scripts/mindlin_test.py'
--- scripts/mindlin_test.py 2010-02-23 14:22:35 +0000
+++ scripts/mindlin_test.py 2010-04-15 16:33:51 +0000
@@ -61,7 +61,7 @@
def myAddPlotData():
i=O.interactions[0,1]
## store some numbers under some labels
- plot.addData(fn=i.phys.normalForce[0],step=O.iter,un=2*s0.shape['radius']-s1.state.pos[0]-s0.state.pos[0],kn=i.phys.kn)
+ plot.addData(fn=i.phys.normalForce[0],step=O.iter,un=2*s0.shape.radius-s1.state.pos[0]+s0.state.pos[0],kn=i.phys.kn)
O.run(50,True);
print "Now calling plot.plot() to show the figure."
