yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2297
committer: Chiara Modenese <chia@engs-018373>
branch nick: trunk
timestamp: Mon 2010-06-21 16:55:15 +0100
message:
  1) Add contact damping into HM law. 2) Add functions getIncidentVel (both for periodic and non-periodic case) and rotateShear to ScGeom. 3) Add py script to demonstrate and test the effect of contact damping.
added:
  scripts/Damping_HM.py
modified:
  pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp*
  pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp*
  pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp*
  pkg/dem/Engine/GlobalEngine/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/DataClass/InteractionGeometry/ScGeom.cpp' (properties changed: -x to +x)
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp	2010-05-26 08:18:36 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp	2010-06-21 15:55:15 +0000
@@ -111,7 +111,7 @@
 		The solution to avoid that is quite simple : use a constant branch vector, like here radius_i*normal.
 		 */
 		// FIXME: For sphere-facet contact this will give an erroneous value of relative velocity...
-		c1x =   radius1*normal; 
+		c1x =   radius1*normal;
 		c2x =  -radius2*normal;
 	} else {
 		// FIXME: It is correct for sphere-sphere and sphere-facet contact
@@ -130,6 +130,71 @@
 	return shearDisplacement;
 }
 
+Vector3r ScGeom::rotateShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt){
+	Vector3r axis;
+	Real angle;
+	// approximated rotations
+		axis = prevNormal.cross(normal); 
+		shearForce -= shearForce.cross(axis);
+		angle = dt*0.5*normal.dot(rbp1->angVel + rbp2->angVel);
+		axis = angle*normal;
+		shearForce -= shearForce.cross(axis);
+	// exact rotations
+	#if 0
+		Quaternionr q;
+		axis					= prevNormal.Cross(normal);
+		angle					= acos(normal.Dot(prevNormal));
+		q.FromAngleAxis(angle,axis);
+		shearForce        = shearForce*q;
+		angle             = dt*0.5*normal.dot(rbp1->angVel+rbp2->angVel);
+		axis					= normal;
+		q.FromAngleAxis(angle,axis);
+		shearForce        = q*shearForce;
+	#endif
+	return shearForce;
+}
+
+Vector3r ScGeom::getIncidentVel(const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting){
+	Vector3r& x = contactPoint;
+	Vector3r c1x, c2x;
+	if(avoidGranularRatcheting){
+		//FIXME : put the long comment on the wiki and keep only a small abstract and link here.
+		/* B.C. Comment : The following definition of c1x and c2x is to avoid "granular ratcheting". This phenomenon has been introduced to me by S. McNamara in a seminar help in Paris, 2004 (GDR MiDi). The concept is also mentionned in many McNamara, Hermann, Lüding, and co-workers papers (see online : "Discrete element methods for the micro-mechanical investigation of granular ratcheting", R. García-Rojo, S. McNamara, H. J. Herrmann, Proceedings ECCOMAS 2004, @ http://www.ica1.uni-stuttgart.de/publications/2004/GMH04/), where it refers to an accumulation of strain in cyclic loadings.
+         
+	        Unfortunately, published papers tend to focus on the "good" ratcheting, i.e. irreversible deformations due to the intrinsic nature of frictional granular materials, while the talk of McNamara in Paris clearly mentioned a possible "bad" ratcheting, purely linked with the formulation of the contact laws in what he called "molecular dynamics" (i.e. Cundall's model, as opposed to "contact dynamics" from Moreau and Jean).
+		
+		Giving a short explanation :
+		The bad ratcheting is best understood considering this small elastic cycle at a contact between two grains : assuming b1 is fixed, impose this displacement to b2 :
+		1. translation "dx" in the normal direction
+		2. rotation "a"
+		3. translation "-dx" (back to initial position)
+		4. rotation "-a" (back to initial orientation)
+		
+		If the branch vector used to define the relative shear in rotation*branch is not constant (typically if it is defined from the vector center->contactPoint like in the "else" below), then the shear displacement at the end of this cycle is not null : rotations a and -a are multiplied by branches of different lengths.
+		It results in a finite contact force at the end of the cycle even though the positions and orientations are unchanged, in total contradiction with the elastic nature of the problem. It could also be seen as an inconsistent energy creation or loss. It is BAD! And given the fact that DEM simulations tend to generate oscillations around equilibrium (damped mass-spring), it can have a significant impact on the evolution of the packings, resulting for instance in slow creep in iterations under constant load.
+		The solution to avoid that is quite simple : use a constant branch vector, like here radius_i*normal.
+		 */
+		// FIXME: For sphere-facet contact this will give an erroneous value of relative velocity...
+		c1x =   radius1*normal;
+		c2x =  -radius2*normal;
+	} else {
+		// FIXME: It is correct for sphere-sphere and sphere-facet contact
+		c1x = (x - rbp1->pos);
+		c2x = (x - rbp2->pos);
+	}
+	Vector3r relativeVelocity = (rbp2->vel+rbp2->angVel.cross(c2x)) - (rbp1->vel+rbp1->angVel.cross(c1x));
+	
+	//update relative velocity for interactions across periods
+	//if (scene->cell->isPeriodic) shearDisplacement+=scene->cell->velGrad*scene->cell->Hsize*cellDist; FIXME : we need pointer to scene. For now, everything will be computed outside this function, which doens't make much sense.
+	relativeVelocity+=shiftVel;
+	return relativeVelocity;
+}
+
+Vector3r ScGeom::getIncidentVel(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){
+	//Just pass null shift to the periodic version
+	return getIncidentVel(rbp1,rbp2,dt,Vector3r::Zero(),avoidGranularRatcheting);
+}
+
 /* keep this for reference; declarations in the header */
 #if 0
 	Vector3r ScGeom::relRotVector() const{

=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp' (properties changed: -x to +x)
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp	2010-05-26 08:18:36 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp	2010-06-21 15:55:15 +0000
@@ -29,6 +29,13 @@
 		//!  Periodic variant. Needs the velocity shift across periods for periodic BCs (else it is safe to pass Vector3r::Zero()). Typically obtained as scene->cell->velGrad*scene->cell->Hsize*cellDist. It would be better to define the shift transparently inside the function, but it needs scene and interaction pointers, which we don't have here.
 		Vector3r rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting=true);
 
+		// Add method which only rotates the shear vector.
+		Vector3r rotateShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt);
+		// Add method which returns the impact velocity (then, inside the contact law, this can be split into shear and normal component). Handle periodicity.
+		Vector3r getIncidentVel(const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting=true);
+		// Implement another version of getIncidentVel which does not handle periodicity.
+		Vector3r getIncidentVel(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
+	
 		YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom,GenericSpheresContact,"Class representing :yref:`geometry<InteractionGeometry>` of two :yref:`spheres<Sphere>` in contact. The contact has 3 DOFs (normal and 2×shear) and uses incremental algorithm for updating shear. (For shear formulated in total displacements and rotations, see :yref:`Dem3DofGeom` and related classes).\n\nWe use symbols $\\vec{x}$, $\\vec{v}$, $\\vec{\\omega}$ respectively for position, linear and angular velocities (all in global coordinates) and $r$ for particles radii; subscripted with 1 or 2 to distinguish 2 spheres in contact. Then we compute unit contact normal\n\n.. math::\n\n\t\\vec{n}=\\frac{\\vec{x}_2-\\vec{x}_1}{||\\vec{x}_2-\\vec{x}_1||}\n\nRelative velocity of spheres is then\n\n.. math::\n\n\t\\vec{v}_{12}=(\\vec{v}_2+\\vec{\\omega}_2\\times(-r_2\\vec{n}))-(\\vec{v}_1+\\vec{\\omega}_1\\times(r_1\\vec{n}))\n\nand its shear component\n\n.. math::\n\n\t\\Delta\\vec{v}_{12}^s=\\vec{v}_{12}-(\\vec{n}\\cdot\\vec{v}_{12})\\vec{n}.\n\nTangential displacement increment over last step then reads\n\n.. math::\n\n\t\\vec{x}_{12}^s=\\Delta t \\vec{v}_{12}^s.",
 		((Vector3r,contactPoint,Vector3r::Zero(),"Reference point of the contact. |ycomp|"))
 		((Vector3r,shear,Vector3r::Zero(),"Total value of the current shear. Update the value using ScGeom::updateShear. |ycomp|"))

=== modified file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp' (properties changed: -x to +x)
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp	2010-06-01 12:04:04 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp	2010-06-21 15:55:15 +0000
@@ -85,32 +85,64 @@
 	/* 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)
+	Real Fn = phys->kno*std::pow(uN,1.5); // normal Force (scalar)
 	phys->normalForce = Fn*scg->normal; // normal Force (vector)
+
+
+	/*** 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){
+		Real mbar = de1->mass*de2->mass / (de1->mass + de2->mass); // equivalent mass
+		Real Cn_crit = 2.*Mathr::Sqrt(mbar*phys->kn); // Critical damping coefficient (normal direction)
+		Real Cs_crit = 2.*Mathr::Sqrt(mbar*phys->ks); // Critical damping coefficient (shear direction)
+		// Note: to compare with the analytical solution you provide cn and cs directly (since here we used a different method to define c_crit)
+		cn = Cn_crit*betan; // Damping normal coefficient
+		cs = Cs_crit*betas; // Damping tangential coefficient
+	}
 
 
 	/*** SHEAR FORCE ***/
-	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)
-	Vector3r& trialFs = phys->shearForce;
-	// define shift to handle periodicity
+	Vector3r& shearElastic = phys->shearElastic; // reference for shearElastic force
+	// Define shift to handle periodicity
 	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 dus =scg->rotateAndGetShear(trialFs, phys->prevNormal, de1, de2, dt, shiftVel, preventGranularRatcheting);
-	//Linear elasticity giving "trial" shear force
-	trialFs -= phys->ks*dus;
-
- 
+	// 1. Rotate shear force
+	shearElastic = scg->rotateShear(shearElastic, phys->prevNormal, de1, de2, dt);
+	// 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, 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 ***/
+	if (useDamping){ // get normal and shear viscous components
+		phys->normalViscous = cn*incidentVn;
+		phys->shearViscous = cs*incidentVs;
+		// add normal viscous component if damping is included
+		phys->normalForce -= phys->normalViscous;
+	}
+
+
 	/*** MOHR-COULOMB LAW ***/
 	Real maxFs = phys->normalForce.squaredNorm()*std::pow(phys->tangensOfFrictionAngle,2);
-	if (trialFs.squaredNorm() > maxFs)
-	{Real ratio = Mathr::Sqrt(maxFs)/trialFs.norm(); trialFs *= ratio;}
+	if (shearElastic.squaredNorm() > maxFs){
+		Real ratio = Mathr::Sqrt(maxFs)/shearElastic.norm(); shearElastic *= ratio; phys->shearForce = shearElastic;}
+	else if (useDamping){ /*add current contact damping if we do not slide and if damping is requested*/ phys->shearForce = shearElastic - phys->shearViscous;}
 
 
 	/*** APPLY FORCES ***/
 	if (!scene->isPeriodic)
-	applyForceAtContactPoint(-phys->normalForce-trialFs , scg->contactPoint , id1, de1->se3.position, id2, de2->se3.position, ncb);
+	applyForceAtContactPoint(-phys->normalForce - phys->shearForce, 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;
+		Vector3r force = -phys->normalForce - phys->shearForce;
 		ncb->forces.addForce(id1,force);
 		ncb->forces.addForce(id2,-force);
 		ncb->forces.addTorque(id1,(scg->radius1-0.5*scg->penetrationDepth)* scg->normal.cross(force));
@@ -119,16 +151,3 @@
 	phys->prevNormal = scg->normal;
 }
 
-
-
-
-
-
-
-
-      
-
-
-	
-	
-

=== modified file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp' (properties changed: -x to +x)
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp	2010-06-01 12:04:04 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp	2010-06-21 15:55:15 +0000
@@ -22,9 +22,13 @@
 class MindlinPhys: public FrictPhys{
 	public:
 	virtual ~MindlinPhys();
-	YADE_CLASS_BASE_DOC_ATTRS_CTOR(MindlinPhys,FrictPhys,"Representation of an interaction of the 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,kso,0.0,"Constant value in the formulation of the tangential stiffness"))
+			((Vector3r,normalViscous,Vector3r::Zero(),"Normal viscous component"))
+			((Vector3r,shearViscous,Vector3r::Zero(),"Normal viscous component"))
+			((Vector3r,shearElastic,Vector3r::Zero(),"Total elastic shear force"))
+			,
 			createIndex());
 	REGISTER_CLASS_INDEX(MindlinPhys,FrictPhys);
 	DECLARE_LOGGER;
@@ -46,11 +50,19 @@
 /******************** Law2_ScGeom_MindlinPhys_Mindlin *********/
 class Law2_ScGeom_MindlinPhys_Mindlin: public LawFunctor{
 	public:
-	virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody);
-	FUNCTOR2D(ScGeom,MindlinPhys);
-	YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_MindlinPhys_Mindlin,LawFunctor,"Constitutive law for the Mindlin's formulation.",
-			((bool,preventGranularRatcheting,true,"bool to avoid granular ratcheting"))
-
+		Real cn, cs;
+		virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody);
+		FUNCTOR2D(ScGeom,MindlinPhys);
+		YADE_CLASS_BASE_DOC_ATTRS_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 contact damping through the definition of the parameters $\\beta_{n}$ and $\\beta_{s}$.",
+			((bool,preventGranularRatcheting,false,"bool to avoid granular ratcheting"))
+			((bool,useDamping,false,"bool to include contact damping"))
+			((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}}$."))
+			,
+			cn = 0.0; cs = 0.0;
+			,
+			.def_readonly("cn",&Law2_ScGeom_MindlinPhys_Mindlin::cn,"Damping normal coefficient.")
+			.def_readonly("cs",&Law2_ScGeom_MindlinPhys_Mindlin::cs,"Damping tangential coefficient.")
 	);
 	DECLARE_LOGGER;
 };

=== added file 'scripts/Damping_HM.py'
--- scripts/Damping_HM.py	1970-01-01 00:00:00 +0000
+++ scripts/Damping_HM.py	2010-06-21 15:55:15 +0000
@@ -0,0 +1,92 @@
+# encoding: utf-8
+# 2010 Chiara Modenese <c.modenese@xxxxxxxxx>
+
+# Script to test the contact damping in HM (both in the normal and shear direction)
+
+from yade import utils
+
+#__________________________________________________________________
+# Geometry sphere and box
+r2=1e-2 # radii [m]
+p2=[r2,0,0] # center positions [m]
+center=[-r2/2.,0,0] # center [m]
+extents=[r2/2.,3*r2,3*r2] # half edge lenght [m]
+
+#__________________________________________________________________
+# Material
+young=600.0e6 # [N/m^2]
+poisson=0.6 
+density=2.60e3 # [kg/m^3]
+frictionAngle=26 # [°]
+
+# Append geometry and material
+O.materials.append(FrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle))
+O.bodies.append(utils.box(center=center,extents=extents,dynamic=False,wire=True)) # body id=0
+O.bodies.append(utils.sphere(p2,r2,dynamic=True,wire=True)) # body id=1
+m_sphere=O.bodies[1].state.mass
+O.bodies[0].state.mass=m_sphere # set the mass of the box the same as the mass of the sphere
+O.bodies[1].state.blockedDOFs=['rx','ry','rz'] # block particles rotations
+
+#__________________________________________________________________
+# list of engines
+O.engines=[
+	ForceResetter(),
+	BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
+	InsertionSortCollider(),
+	InteractionDispatchers(
+		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
+		[Ip2_FrictMat_FrictMat_MindlinPhys(label='damping')],
+		[Law2_ScGeom_MindlinPhys_Mindlin(label='contactLaw')]
+	),
+	ForceEngine(force=(-20,0,0),subscribedBodies=[1],label='force'),	
+	NewtonIntegrator(damping=0.0),
+	PeriodicPythonRunner(iterPeriod=1,command='myAddPlotData()'),
+]
+
+#__________________________________________________________________
+# define contact damping coefficients
+contactLaw.betan=0.5 # normal direction
+contactLaw.betas=0.0 # shear direction
+contactLaw.useDamping=True
+
+#__________________________________________________________________
+# time step
+O.dt=.2*utils.PWaveTimeStep()
+O.saveTmp('init')
+
+#__________________________________________________________________
+from yade import qt
+qt.View()
+qt.Controller()
+
+#__________________________________________________________________
+# plot some results
+from math import *
+from yade import plot
+
+plot.labels=dict(Fn='$Normal\,force$',un='$Overlapping$',time='$Time$',time_='$Time$',Fs='$Shear\,force$',t='$Time$')
+plot.plots={'time':('un',),'time_':('Fn',),'t':('Fs')}
+
+def myAddPlotData():
+	i=O.interactions[0,1]
+	plot.addData(Fn=i.phys.normalForce[0],Fs=i.phys.shearForce[1],un=i.geom.penetrationDepth,time=O.time,time_=O.time,t=O.time)
+
+# We will have:
+# 1) data in graphs (if you call plot.plot())
+# 2) data in file (if you call plot.saveGnuplot('/home/chia/Desktop/Output/out')
+# 3) data in memory as plot.data['un'], plot.data['fn'], etc. under the labels they were saved
+
+O.run(1000,True)
+
+# Now test the shear direction, having previously obtained a constant normal overlapping 
+contactLaw.betas=0.5
+force.force=(-20,5,0) # assign a force also in the shear direction
+O.run(1000,True)
+plot.plot()
+
+# NOTE about the results:
+# In the graphs, you will see the behaviour both of the normal and shear force. As expected from the damped solution, their values oscillate around the equilibrium position, which is eventually reached after few iterations. 
+
+
+
+