yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3984
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
timestamp: Mon 2014-05-26 18:49:23 +0200
message:
  CohesiveFrictional contacts: make them elastic only if strength<0 (previous condition was <=0) + Lingran's elastic energy of contact moments
modified:
  pkg/dem/CohFrictMat.hpp
  pkg/dem/CohFrictPhys.hpp
  pkg/dem/CohesiveFrictionalContactLaw.cpp
  pkg/dem/CohesiveFrictionalContactLaw.hpp
  pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== modified file 'pkg/dem/CohFrictMat.hpp'
--- pkg/dem/CohFrictMat.hpp	2014-05-19 15:00:11 +0000
+++ pkg/dem/CohFrictMat.hpp	2014-05-26 16:49:23 +0000
@@ -23,7 +23,8 @@
 		((bool,isCohesive,true,,""))
 		((Real,alphaKr,2.0,,"Dimensionless rolling stiffness."))
 		((Real,alphaKtw,2.0,,"Dimensionless twist stiffness."))
-		((Real,etaRoll,-1.,,"Dimensionless rolling (aka 'bending') strength. If negative, rolling will be elastic."))	
+		((Real,etaRoll,-1.,,"Dimensionless rolling (aka 'bending') strength. If negative, rolling moment will be elastic."))
+		((Real,etaTwist,-1.,,"Dimensionless twisting strength. If negative, twist moment will be elastic."))
 		((Real,normalCohesion,0,,"Tensile strength, homogeneous to a pressure"))
 		((Real,shearCohesion,0,,"Shear strength, homogeneous to a pressure"))
 		((bool,momentRotationLaw,false,,"Use bending/twisting moment at contact. The contact will have moments only if both bodies have this flag true. See :yref:`CohFrictPhys::cohesionDisablesFriction` for details."))

=== modified file 'pkg/dem/CohFrictPhys.hpp'
--- pkg/dem/CohFrictPhys.hpp	2014-01-21 15:08:51 +0000
+++ pkg/dem/CohFrictPhys.hpp	2014-05-26 16:49:23 +0000
@@ -23,8 +23,8 @@
 		((bool,fragile,true,,"do cohesion disapear when contact strength is exceeded?"))
 		((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,maxRollPl,0.0,,"Coefficient of rolling friction (negative means elastic)."))
+		((Real,maxTwistPl,0.0,,"Coefficient of twisting friction (negative means elastic)."))
 		((Real,normalAdhesion,0,,"tensile strength"))
 		((Real,shearAdhesion,0,,"cohesive part of the shear strength (a frictional term might be added depending on :yref:`CohFrictPhys::cohesionDisablesFriction`)"))
 		((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	2013-03-06 17:30:45 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.cpp	2014-05-26 16:49:23 +0000
@@ -16,6 +16,9 @@
 YADE_PLUGIN((CohesiveFrictionalContactLaw)(Law2_ScGeom6D_CohFrictPhys_CohesionMoment));
 CREATE_LOGGER(Law2_ScGeom6D_CohFrictPhys_CohesionMoment);
 
+Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::getPlasticDissipation() {return (Real) plasticDissipation;}
+void Law2_ScGeom6D_CohFrictPhys_CohesionMoment::initPlasticDissipation(Real initVal) {plasticDissipation.reset(); plasticDissipation+=initVal;}
+
 Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::normElastEnergy()
 {
 	Real normEnergy=0;
@@ -41,6 +44,50 @@
 	return shearEnergy;
 }
 
+Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::bendingElastEnergy()
+{
+	Real bendingEnergy=0;
+	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+		if(!I->isReal()) continue;
+		CohFrictPhys* phys = YADE_CAST<CohFrictPhys*>(I->phys.get());
+		if (phys) {
+			bendingEnergy += 0.5*(phys->moment_bending.squaredNorm()/phys->kr);
+		}
+	}
+	return bendingEnergy;
+}
+
+Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::twistElastEnergy()
+{
+	Real twistEnergy=0;
+	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+		if(!I->isReal()) continue;
+		CohFrictPhys* phys = YADE_CAST<CohFrictPhys*>(I->phys.get());
+		if (phys) {
+			twistEnergy += 0.5*(phys->moment_twist.squaredNorm()/phys->ktw);
+		}
+	}
+	return twistEnergy;
+}
+
+Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::totalElastEnergy()
+{
+	Real totalEnergy=0;
+	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+		if(!I->isReal()) continue;
+		CohFrictPhys* phys = YADE_CAST<CohFrictPhys*>(I->phys.get());
+		if (phys) {
+			totalEnergy += 0.5*(phys->normalForce.squaredNorm()/phys->kn);
+			totalEnergy += 0.5*(phys->shearForce.squaredNorm()/phys->ks);
+			totalEnergy += 0.5*(phys->moment_bending.squaredNorm()/phys->kr);
+			totalEnergy += 0.5*(phys->moment_twist.squaredNorm()/phys->ktw);
+		}
+	}
+	return totalEnergy;
+}
+
+
+
 void CohesiveFrictionalContactLaw::action()
 {
 	if(!functor) functor=shared_ptr<Law2_ScGeom6D_CohFrictPhys_CohesionMoment>(new Law2_ScGeom6D_CohFrictPhys_CohesionMoment);
@@ -104,17 +151,12 @@
 			Vector3r trialForce=shearForce;
 			shearForce *= maxFs;
 			if (scene->trackEnergy){
-				Real dissip=((1/phys->ks)*(trialForce-shearForce))/*plastic disp*/ .dot(shearForce)/*active force*/;
-				if(dissip>0) scene->energy->add(dissip,"plastDissip",plastDissipIx,/*reset*/false);}
+				Real sheardissip=((1/phys->ks)*(trialForce-shearForce))/*plastic disp*/ .dot(shearForce)/*active force*/;
+				if(sheardissip>0) scene->energy->add(sheardissip,"shearDissip",shearDissipIx,/*reset*/false);}
 			if (Fn<0)  phys->normalForce = Vector3r::Zero();//Vector3r::Zero()
 		}
 		//Apply the force
 		applyForceAtContactPoint(-phys->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position + (scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero()));
-// 		Vector3r force = -phys->normalForce-shearForce;
-// 		scene->forces.addForce(id1,force);
-// 		scene->forces.addForce(id2,-force);
-// 		scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
-// 		scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
 
 		/// Moment law  ///
 		if (phys->momentRotationLaw && (!phys->cohesionBroken || always_use_moment_law)) {
@@ -152,22 +194,28 @@
 			/// Plasticity ///
 			// limit rolling moment to the plastic value, if required
 			Real RollMax = phys->maxRollPl*phys->normalForce.norm();
-			if (RollMax>0.){ // do we want to apply plasticity?
+			if (RollMax>=0.){ // do we want to apply plasticity?
 				if (!useIncrementalForm) LOG_WARN("If :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::useIncrementalForm` is false, then plasticity will not be applied correctly (the total formulation would not reproduce irreversibility).");
 				Real scalarRoll = phys->moment_bending.norm();
 				if (scalarRoll>RollMax){ // fix maximum rolling moment
 					Real ratio = RollMax/scalarRoll;
 					phys->moment_bending *= ratio;
-				}	
+					if (scene->trackEnergy){
+						Real bendingdissip=((1/phys->kr)*(scalarRoll-RollMax)*RollMax)/*active force*/;
+						if(bendingdissip>0) scene->energy->add(bendingdissip,"bendingDissip",bendingDissipIx,/*reset*/false);}
+				}
 			}
 			// limit twisting moment to the plastic value, if required
-			Real TwistMax = phys->maxTwistMoment.norm();
-			if (TwistMax>0.){ // do we want to apply plasticity?
+			Real TwistMax = phys->maxTwistPl*phys->normalForce.norm();
+			if (TwistMax>=0.){ // do we want to apply plasticity?
 				if (!useIncrementalForm) LOG_WARN("If :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::useIncrementalForm` is false, then plasticity will not be applied correctly (the total formulation would not reproduce irreversibility).");
 				Real scalarTwist= phys->moment_twist.norm();
 				if (scalarTwist>TwistMax){ // fix maximum rolling moment
 					Real ratio = TwistMax/scalarTwist;
 					phys->moment_twist *= ratio;
+					if (scene->trackEnergy){
+						Real twistdissip=((1/phys->ktw)*(scalarTwist-TwistMax)*TwistMax)/*active force*/;
+						if(twistdissip>0) scene->energy->add(twistdissip,"twistDissip",twistDissipIx,/*reset*/false);}
 				}	
 			}
 			// Apply moments now

=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/CohesiveFrictionalContactLaw.hpp	2014-01-21 15:06:34 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.hpp	2014-05-26 16:49:23 +0000
@@ -17,20 +17,32 @@
 
 class Law2_ScGeom6D_CohFrictPhys_CohesionMoment: public LawFunctor{
 	public:
+		OpenMPAccumulator<Real> plasticDissipation;
 		Real normElastEnergy();
 		Real shearElastEnergy();
+		Real bendingElastEnergy();
+		Real twistElastEnergy();
+		Real totalElastEnergy();
+		Real getPlasticDissipation();
+		void initPlasticDissipation(Real initVal=0);
 	virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
 	YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_CohFrictPhys_CohesionMoment,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_s$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False (in course of implementation), the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$ (the two values are the same currently), so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies. The maximum value of moments can be defined and takes the form of rolling friction. Cohesive -type moment may also be included in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=-F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=-M_t/\\nu_t$.",
 		((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
 		((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,traceEnergy,false,,"Define the total energy dissipated in plastic slips at all contacts. This will trace only plastic energy in this law, see O.trackEnergy for a more complete energies tracing"))
 		((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."))
-		((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)"))
+		((int,shearDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for shear dissipation (with O.trackEnergy)"))
+		((int,bendingDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for bending dissipation (with O.trackEnergy)"))
+		((int,twistDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for twist dissipation (with O.trackEnergy)"))
 		((Real,creep_viscosity,1,,"creep viscosity [Pa.s/m]. probably should be moved to Ip2_CohFrictMat_CohFrictMat_CohFrictPhys."))
 		,,
 		.def("normElastEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::normElastEnergy,"Compute normal elastic energy.")
 		.def("shearElastEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::shearElastEnergy,"Compute shear elastic energy.")
+		.def("bendingElastEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::bendingElastEnergy,"Compute bending elastic energy.")
+		.def("twistElastEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::twistElastEnergy,"Compute twist elastic energy.")
+		.def("elasticEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::totalElastEnergy,"Compute total elastic energy.")
 	);
 	FUNCTOR2D(ScGeom6D,CohFrictPhys);
 	DECLARE_LOGGER;

=== modified file 'pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp'
--- pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp	2014-01-24 17:15:20 +0000
+++ pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp	2014-05-26 16:49:23 +0000
@@ -51,9 +51,6 @@
 			if (Va && Vb) Ks = 2.0*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Vb);//harmonic average of two stiffnesses with ks=V*kn for each sphere
 			else Ks=0;
 
-			// Jean-Patrick Plassiard, Noura Belhaine, Frederic
-			// Victor Donze, "Calibration procedure for spherical
-			// discrete elements using a local moemnt law".
 			contactPhysics->kr = Da*Db*Ks*AlphaKr;
 			contactPhysics->ktw = Da*Db*Ks*AlphaKtw;
 			contactPhysics->tangensOfFrictionAngle		= std::tan(std::min(fa,fb));
@@ -69,8 +66,8 @@
 			contactPhysics->ks = Ks;
 
 			contactPhysics->maxRollPl = min(sdec1->etaRoll*Da,sdec2->etaRoll*Db);
+			contactPhysics->maxTwistPl = min(sdec1->etaTwist*Da,sdec2->etaTwist*Db);
 			contactPhysics->momentRotationLaw=(sdec1->momentRotationLaw && sdec2->momentRotationLaw);
-			//contactPhysics->elasticRollingLimit = elasticRollingLimit;
 		}
 		else {// !isNew, but if setCohesionNow, all contacts are initialized like if they were newly created
 			CohFrictPhys* contactPhysics = YADE_CAST<CohFrictPhys*>(interaction->phys.get());