yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2681
committer: Chiara Modenese <c.modenese@xxxxxxxxx>
branch nick: yade
timestamp: Wed 2011-01-26 18:19:52 +0000
message:
  - Add the possibility to calculate bending and twisting moments using the incremental formulation in both CFLaw and HMLaw;
  - Add maximum elastic torsional moment in CFLaw (it should be computed using normal force, I will see how to compute it somehow);
  - Fix mistake in mindlin.py example script.
modified:
  pkg/dem/CohFrictPhys.hpp
  pkg/dem/CohesiveFrictionalContactLaw.cpp
  pkg/dem/CohesiveFrictionalContactLaw.hpp
  pkg/dem/HertzMindlin.cpp
  pkg/dem/HertzMindlin.hpp
  scripts/test/mindlin.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/CohFrictPhys.hpp'
--- pkg/dem/CohFrictPhys.hpp	2010-12-02 19:27:56 +0000
+++ pkg/dem/CohFrictPhys.hpp	2011-01-26 18:19:52 +0000
@@ -24,6 +24,7 @@
 		((Real,kr,0,,"rotational stiffness [N.m/rad]"))
 		((Real,ktw,0,,"twist stiffness [N.m/rad]"))
 		((Real,maxRollPl,0.0,,"Coefficient to determine the maximum plastic rolling moment."))
+		((Vector3r,maxTwistMoment,Vector3r::Zero(),,"Maximum elastic value for the twisting moment (if zero, plasticity will not be applied). In CohFrictMat a parameter should be added to decide what value should be attributed to this threshold value."))
 		((Real,normalAdhesion,0,,"tensile strength"))
 		((Real,shearAdhesion,0,,"cohesive part of the shear strength (a frictional term might be added depending on :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::always_use_moment_law`)"))
 		((Real,unp,0,,"plastic normal displacement, only used for tensile behaviour and if :yref:`CohFrictPhys::fragile`=false."))

=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/CohesiveFrictionalContactLaw.cpp	2010-12-08 10:43:28 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.cpp	2011-01-26 18:19:52 +0000
@@ -83,29 +83,61 @@
 		}
 		applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
 
-		/// Moment law        ///
+		/// Moment law  ///
 		if (currentContactPhysics->momentRotationLaw && (!currentContactPhysics->cohesionBroken || always_use_moment_law)) {
-			if (twist_creep) {
-				Real viscosity_twist = creep_viscosity * std::pow((2 * std::min(currentContactGeometry->radius1,currentContactGeometry->radius2)),2) / 16.0;
-				Real angle_twist_creeped = currentContactGeometry->getTwist() * (1 - dt/viscosity_twist);
-				Quaternionr q_twist(AngleAxisr(currentContactGeometry->getTwist(),currentContactGeometry->normal));
-				Quaternionr q_twist_creeped(AngleAxisr(angle_twist_creeped,currentContactGeometry->normal));
-				Quaternionr q_twist_delta(q_twist_creeped * q_twist.conjugate());
-				currentContactGeometry->twistCreep = currentContactGeometry->twistCreep * q_twist_delta;
+			if (!useIncrementalForm){
+				if (twist_creep) {
+					Real viscosity_twist = creep_viscosity * std::pow((2 * std::min(currentContactGeometry->radius1,currentContactGeometry->radius2)),2) / 16.0;
+					Real angle_twist_creeped = currentContactGeometry->getTwist() * (1 - dt/viscosity_twist);
+					Quaternionr q_twist(AngleAxisr(currentContactGeometry->getTwist(),currentContactGeometry->normal));
+					Quaternionr q_twist_creeped(AngleAxisr(angle_twist_creeped,currentContactGeometry->normal));
+					Quaternionr q_twist_delta(q_twist_creeped * q_twist.conjugate());
+					currentContactGeometry->twistCreep = currentContactGeometry->twistCreep * q_twist_delta;
+				}
+				currentContactPhysics->moment_twist = (currentContactGeometry->getTwist()*currentContactPhysics->ktw)*currentContactGeometry->normal;
+				currentContactPhysics->moment_bending = currentContactGeometry->getBending() * currentContactPhysics->kr;
+			}	
+			else{ // Use incremental formulation to compute moment_twis and moment_bending (no twist_creep is applied)
+				if (twist_creep) throw std::invalid_argument("Law2_ScGeom6D_CohFrictPhys_CohesionMoment: no twis creep is included if the incremental form for the rotations is used.");
+				Vector3r relAngVel = currentContactGeometry->getRelAngVel(de1,de2,dt);
+				// *** Bending ***//
+				Vector3r relAngVelBend = relAngVel - currentContactGeometry->normal.dot(relAngVel)*currentContactGeometry->normal; // keep only the bending part 
+				Vector3r relRotBend = relAngVelBend*dt; // relative rotation due to rolling behaviour	
+				// incremental formulation for the bending moment (as for the shear part)
+				Vector3r& momentBend = currentContactPhysics->moment_bending;
+				momentBend = currentContactGeometry->rotate(momentBend); // rotate moment vector (updated)
+				momentBend = momentBend-currentContactPhysics->kr*relRotBend;
+				// ----------------------------------------------------------------------------------------
+				// *** Torsion ***//
+				Vector3r relAngVelTwist = currentContactGeometry->normal.dot(relAngVel)*currentContactGeometry->normal;
+				Vector3r relRotTwist = relAngVelTwist*dt; // component of relative rotation along n  FIXME: sign?
+				// incremental formulation for the torsional moment
+				Vector3r& momentTwist = currentContactPhysics->moment_twist;
+				momentTwist = currentContactGeometry->rotate(momentTwist); // rotate moment vector (updated)
+				momentTwist = momentTwist-currentContactPhysics->ktw*relRotTwist; // FIXME: sign?
 			}
-			currentContactPhysics->moment_twist = (currentContactGeometry->getTwist()*currentContactPhysics->ktw)*currentContactGeometry->normal;
-			currentContactPhysics->moment_bending = currentContactGeometry->getBending() * currentContactPhysics->kr;
-			
+			/// Plasticity ///
 			// limit rolling moment to the plastic value, if required
 			Real RollMax = currentContactPhysics->maxRollPl*currentContactPhysics->normalForce.norm();
 			if (RollMax>0.){ // do we want to apply plasticity?
+				LOG_WARN("If :yref:`CohesiveFrictionalContactLaw::useIncrementalForm` is false, then plasticity would not be applied correctly (the total formulation would not reproduce irreversibility).");
 				Real scalarRoll = currentContactPhysics->moment_bending.norm();		
 				if (scalarRoll>RollMax){ // fix maximum rolling moment
 					Real ratio = RollMax/scalarRoll;
-					currentContactPhysics->moment_bending *= ratio;		
-				}	
-			}
-			
+					currentContactPhysics->moment_bending *= ratio;	
+				}	
+			}
+			// limit twisting moment to the plastic value, if required
+			Real TwistMax = currentContactPhysics->maxTwistMoment.norm();
+			if (TwistMax>0.){ // do we want to apply plasticity?
+				LOG_WARN("If :yref:`CohesiveFrictionalContactLaw::useIncrementalForm` is false, then plasticity would not be applied correctly (the total formulation would not reproduce irreversibility).");
+				Real scalarTwist= currentContactPhysics->moment_twist.norm();		
+				if (scalarTwist>TwistMax){ // fix maximum rolling moment
+					Real ratio = TwistMax/scalarTwist;
+					currentContactPhysics->moment_twist *= ratio;	
+				}	
+			}
+			// Apply moments now
 			Vector3r moment = currentContactPhysics->moment_twist + currentContactPhysics->moment_bending;
 			scene->forces.addTorque(id1,-moment);
 			scene->forces.addTorque(id2, moment);			

=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/CohesiveFrictionalContactLaw.hpp	2011-01-09 16:34:50 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.hpp	2011-01-26 18:19:52 +0000
@@ -23,10 +23,13 @@
 		((bool,always_use_moment_law,false,,"If true, use bending/twisting moments at all contacts. If false, compute moments only for cohesive contacts."))
 		((bool,shear_creep,false,,"activate creep on the shear force, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
 		((bool,twist_creep,false,,"activate creep on the twisting moment, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
+		((bool,useIncrementalForm,false,,"use the incremental formulation to compute bending and twisting moments. Creep on the twisting moment is not included in such a case."))
 		((Real,creep_viscosity,1,,"creep viscosity [Pa.s/m]. probably should be moved to Ip2_CohFrictMat_CohFrictMat_CohFrictPhys..."))
 	);
 	FUNCTOR2D(ScGeom6D,CohFrictPhys);
 	DECLARE_LOGGER;
+	private:
+	Real currentContactPhysics;
 };
 REGISTER_SERIALIZABLE(Law2_ScGeom6D_CohFrictPhys_CohesionMoment);
 

=== modified file 'pkg/dem/HertzMindlin.cpp'
--- pkg/dem/HertzMindlin.cpp	2011-01-21 19:55:58 +0000
+++ pkg/dem/HertzMindlin.cpp	2011-01-26 18:19:52 +0000
@@ -76,6 +76,7 @@
 	mindlinPhys->adhesionForce = Adhesion;
 	
 	mindlinPhys->kr = krot;
+	mindlinPhys->ktw = ktwist;
 	mindlinPhys->maxBendPl = eta*Rmean; // does this make sense? why do we take Rmean?
 
 	/* compute viscous coefficients */
@@ -411,19 +412,28 @@
 	}
 	
 	/********************************************/
-	/* MOMENT CONTACT LAW (only bending moment) */
+	/* MOMENT CONTACT LAW */
 	/********************************************/
 	if (includeMoment){
+		// *** Bending ***//
 		// new code to compute relative particle rotation (similar to the way the shear is computed)
 		// use scg function to compute relAngVel
 		Vector3r relAngVel = scg->getRelAngVel(de1,de2,dt);
 		//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	
+		Vector3r relAngVelBend = relAngVel - scg->normal.dot(relAngVel)*scg->normal; // keep only the bending part 
+		Vector3r relRot = relAngVelBend*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?
+		momentBend = momentBend-phys->kr*relRot; // add incremental rolling to the rolling vector
+		// ----------------------------------------------------------------------------------------
+		// *** Torsion ***//
+		Vector3r relAngVelTwist = scg->normal.dot(relAngVel)*scg->normal;
+		Vector3r relRotTwist = relAngVelTwist*dt; // component of relative rotation along n
+		// incremental formulation for the torsional moment
+		Vector3r& momentTwist = phys->momentTwist;
+		momentTwist = scg->rotate(momentTwist); // rotate moment vector (updated)
+		momentTwist = momentTwist-phys->ktw*relRotTwist;
 
 #if 0
 	// code to compute the relative particle rotation
@@ -455,7 +465,7 @@
 		momentBend = momentBend+phys->kr*phys->dThetaR; // add incremental rolling to the rolling vector FIXME: is the sign correct?
 #endif
 
-		// check plasticity condition
+		// check plasticity condition (only bending part for the moment)
 		Real MomentMax = phys->maxBendPl*phys->normalForce.norm();
 		Real scalarMoment = phys->momentBend.norm();
 		if (MomentMax>0){
@@ -466,8 +476,9 @@
 			 }
 		}
 		// apply moments
-		scene->forces.addTorque(id1,-phys->momentBend); 
-		scene->forces.addTorque(id2,phys->momentBend);
+		Vector3r moment = phys->momentTwist+phys->momentBend;
+		scene->forces.addTorque(id1,-moment); 
+		scene->forces.addTorque(id2,moment);
 	}
 
 	// update variables

=== modified file 'pkg/dem/HertzMindlin.hpp'
--- pkg/dem/HertzMindlin.hpp	2011-01-13 16:46:20 +0000
+++ pkg/dem/HertzMindlin.hpp	2011-01-26 18:19:52 +0000
@@ -32,6 +32,7 @@
 			((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,ktw,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"))
@@ -42,6 +43,7 @@
 
 			//((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,momentTwist,Vector3r::Zero(),,"Artificial twisting moment (no plastic condition can be applied at the moment)"))
 			//((Vector3r,dThetaR,Vector3r::Zero(),,"Incremental rolling vector"))
 
 			((Real,radius,NaN,,"Contact radius (only computed with :yref:`Law2_ScGeom_MindlinPhys_Mindlin::calcEnergy`)"))
@@ -78,6 +80,7 @@
 			((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 determine the plastic bending moment"))
 			((Real,krot,0.0,,"Rotational stiffness for moment contact law"))
+			((Real,ktwist,0.0,,"Torsional 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$."))

=== modified file 'scripts/test/mindlin.py'
--- scripts/test/mindlin.py	2010-11-07 08:55:43 +0000
+++ scripts/test/mindlin.py	2011-01-26 18:19:52 +0000
@@ -12,7 +12,7 @@
 		[Ip2_FrictMat_FrictMat_MindlinPhys()],
 		[Law2_ScGeom_MindlinPhys_Mindlin()]
 	),
-	GravityEngine(gravity=(-100,0,0)),
+	GravityEngine(gravity=(10,0,0)),
 	NewtonIntegrator(damping=0.0),
 	###
 	### NOTEÂ this extra engine:
@@ -28,8 +28,8 @@
 
 ## create two spheres (one will be fixed) and append them
 from yade import utils
-s0=utils.sphere([0,0,0],1,color=[0,1,0],dynamic=False,wire=True,material='Friction')
-s1=utils.sphere([2,0,0],1,color=[0,2,0],dynamic=True,wire=True,material='Friction')
+s0=utils.sphere([0,0,0],1,color=[0,1,0],fixed=False,wire=True,material='Friction')
+s1=utils.sphere([2,0,0],1,color=[0,2,0],fixed=True,wire=True,material='Friction')
 O.bodies.append(s0)
 O.bodies.append(s1)
 
@@ -59,8 +59,7 @@
 	## 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)	
 
-O.run(250,True);
-print "Now calling plot.plot() to show the figure."
+O.run(100,True); plot.plot()
 
 ## We will have:
 ## 1) data in graphs (if you call plot.plot())