yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3695
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
timestamp: Wed 2015-07-01 20:36:34 +0200
message:
  account for the rotational stiffness of interactions in GlobalStiffnessTimeStepper
modified:
  pkg/dem/CohesiveFrictionalContactLaw.cpp
  pkg/dem/CohesiveFrictionalContactLaw.hpp
  pkg/dem/GlobalStiffnessTimeStepper.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/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/CohesiveFrictionalContactLaw.cpp	2015-02-12 11:14:24 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.cpp	2015-07-01 18:36:34 +0000
@@ -17,6 +17,9 @@
 Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::getPlasticDissipation() {return (Real) plasticDissipation;}
 void Law2_ScGeom6D_CohFrictPhys_CohesionMoment::initPlasticDissipation(Real initVal) {plasticDissipation.reset(); plasticDissipation+=initVal;}
 
+Vector3r CohFrictPhys::getRotStiffness() {return Vector3r(ktw,kr,kr);}
+
+
 Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::normElastEnergy()
 {
 	Real normEnergy=0;
@@ -256,8 +259,9 @@
 			Real Kn = 2.0*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
 
 			// harmonic average of alphas parameters
-			Real AlphaKr = 2.0*sdec1->alphaKr*sdec2->alphaKr/(sdec1->alphaKr+sdec2->alphaKr);
-			Real AlphaKtw;
+			Real AlphaKr, AlphaKtw;
+			if (sdec1->alphaKr && sdec2->alphaKr) AlphaKr = 2.0*sdec1->alphaKr*sdec2->alphaKr/(sdec1->alphaKr+sdec2->alphaKr);
+			else AlphaKr = 0;
 			if (sdec1->alphaKtw && sdec2->alphaKtw) AlphaKtw = 2.0*sdec1->alphaKtw*sdec2->alphaKtw/(sdec1->alphaKtw+sdec2->alphaKtw);
 			else AlphaKtw=0;
 

=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/CohesiveFrictionalContactLaw.hpp	2015-04-30 14:43:53 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.hpp	2015-07-01 18:36:34 +0000
@@ -22,7 +22,6 @@
 {
 	public :
 		virtual ~CohFrictMat () {};
-
 /// Serialization
 	YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohFrictMat,FrictMat,"",
 		((bool,isCohesive,true,,""))
@@ -49,6 +48,7 @@
 	public :
 		virtual ~CohFrictPhys() {};
 		void SetBreakingState() {cohesionBroken = true; normalAdhesion = 0; shearAdhesion = 0;};
+		virtual Vector3r getRotStiffness();
 
 	YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohFrictPhys,FrictPhys,"",
 		((bool,cohesionDisablesFriction,false,,"is shear strength the sum of friction and adhesion or only adhesion?"))

=== modified file 'pkg/dem/GlobalStiffnessTimeStepper.cpp'
--- pkg/dem/GlobalStiffnessTimeStepper.cpp	2014-10-15 06:44:01 +0000
+++ pkg/dem/GlobalStiffnessTimeStepper.cpp	2015-07-01 18:36:34 +0000
@@ -177,13 +177,16 @@
 				std::pow(normal.x(),2)+std::pow(normal.z(),2),
 				std::pow(normal.x(),2)+std::pow(normal.y(),2));
 		diag_Rstiffness *= ks;
-
-
-		//NOTE : contact laws with moments would be handled correctly by summing directly bending+twisting stiffness to diag_Rstiffness. The fact that there is no problem currently with e.g. cohesiveFrict law is probably because final computed dt is constrained by translational motion, not rotations.
+		
+		// If contact moments are present, add the diagonal of (n⊗n*k_twist + (I-n⊗n)*k_roll = (k_twist-k_roll)*n⊗n + I*k_roll ) :
+		Vector3r kr = (static_cast<NormShearPhys *> (contact->phys.get()))->getRotStiffness(); //get the vector (k_twist,k_roll,k_roll)
+		Vector3r nn (std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));//n⊗n 
+		Vector3r diag_Mstiffness= (kr[0]-kr[1])*nn + Vector3r(1,1,1)*kr[1];
+
 		stiffnesses [contact->getId1()]+=diag_stiffness;
-		Rstiffnesses[contact->getId1()]+=diag_Rstiffness*pow(radius1,2);
+		Rstiffnesses[contact->getId1()]+=diag_Rstiffness*pow(radius1,2)+ diag_Mstiffness;
 		stiffnesses [contact->getId2()]+=diag_stiffness;
-		Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);	
+		Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2)+ diag_Mstiffness;
 
 		//Same for the Viscous part, if required
 		if (viscEl == true){