yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

Merge authors:
  help
------------------------------------------------------------
revno: 2596 [merge]
committer: Chiara Modenese <c.modenese@xxxxxxxxx>
branch nick: yade
timestamp: Wed 2010-12-08 11:29:11 +0000
message:
  1. Move alphas parameters to FrictMat as discussed with Bruno (harmonic-average of them is introduced in Ip2 functor).
   2. Forgot one inequality when plasticity is applied (maybe I do not like tests so much :-) Sorry, anyway).
modified:
  pkg/dem/CohFrictMat.hpp
  pkg/dem/CohesiveFrictionalContactLaw.cpp
  pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.cpp
  pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.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/CohFrictMat.hpp'
--- pkg/dem/CohFrictMat.hpp	2010-12-06 19:24:20 +0000
+++ pkg/dem/CohFrictMat.hpp	2010-12-08 10:43:28 +0000
@@ -21,6 +21,9 @@
 /// Serialization
 	YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohFrictMat,FrictMat,"",
 		((bool,isCohesive,true,,""))
+		((Real,alphaKr,2.0,,"Dimensionless coefficient used for the rolling stiffness."))
+		((Real,alphaKtw,2.0,,"Dimensionless coefficient used for the twist stiffness."))
+		((Real,etaRoll,-1.,,"Dimensionless coefficient used to calculate the plastic rolling moment (if negative, plasticity will not be applied)."))	
 		((Real,normalCohesion,0,,""))
 		((Real,shearCohesion,0,,""))
 		((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/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/CohesiveFrictionalContactLaw.cpp	2010-12-02 19:27:56 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.cpp	2010-12-08 10:43:28 +0000
@@ -98,10 +98,12 @@
 			
 			// limit rolling moment to the plastic value, if required
 			Real RollMax = currentContactPhysics->maxRollPl*currentContactPhysics->normalForce.norm();
-			if (RollMax>0.){ // apply plasticity
+			if (RollMax>0.){ // do we want to apply plasticity?
 				Real scalarRoll = currentContactPhysics->moment_bending.norm();		
-				Real ratio = RollMax/scalarRoll;
-				currentContactPhysics->moment_bending *= ratio;		
+				if (scalarRoll>RollMax){ // fix maximum rolling moment
+					Real ratio = RollMax/scalarRoll;
+					currentContactPhysics->moment_bending *= ratio;		
+				}	
 			}
 			
 			Vector3r moment = currentContactPhysics->moment_twist + currentContactPhysics->moment_bending;

=== modified file 'pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.cpp'
--- pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.cpp	2010-12-03 12:00:34 +0000
+++ pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.cpp	2010-12-08 10:43:28 +0000
@@ -39,11 +39,15 @@
 			Real Vb 	= sdec2->poisson;
 			Real Da 	= geom->radius1;
 			Real Db 	= geom->radius2;
-			Real Dmean  = (Da+Db)/2.; // mean radius
 			Real fa 	= sdec1->frictionAngle;
 			Real fb 	= sdec2->frictionAngle;
 			Real Kn = 2.0*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
 
+			// harmonic average of alphas and etaRoll parameters
+			Real AlphaKr = 2.0*sdec1->alphaKr*sdec2->alphaKr/(sdec1->alphaKr+sdec2->alphaKr);
+			Real AlphaKtw = 2.0*sdec1->alphaKtw*sdec2->alphaKtw/(sdec1->alphaKtw+sdec2->alphaKtw);
+			Real EtaRoll = 2.0*sdec1->etaRoll*sdec2->etaRoll/(sdec1->etaRoll+sdec2->etaRoll);
+
 			Real Ks;
 			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;
@@ -51,13 +55,8 @@
 			// Jean-Patrick Plassiard, Noura Belhaine, Frederic
 			// Victor Donze, "Calibration procedure for spherical
 			// discrete elements using a local moemnt law".
-			if (!useMeanRad){
-				contactPhysics->kr = Da*Db*Ks*alphaKr;
-				contactPhysics->ktw = Da*Db*Ks*alphaKtw;
-			} else {
-			// Chiara Modese's variant
-				contactPhysics->kr = pow(Dmean,2)*Ks*alphaKr;
-				contactPhysics->ktw = pow(Dmean,2)*Ks*alphaKtw;}
+			contactPhysics->kr = Da*Db*Ks*AlphaKr;
+			contactPhysics->ktw = Da*Db*Ks*AlphaKtw;
 
 			contactPhysics->frictionAngle			= std::min(fa,fb);
 			contactPhysics->tangensOfFrictionAngle		= std::tan(contactPhysics->frictionAngle);
@@ -72,7 +71,7 @@
 			contactPhysics->kn = Kn;
 			contactPhysics->ks = Ks;
 
-			contactPhysics->maxRollPl = etaRoll*Dmean;
+			contactPhysics->maxRollPl = EtaRoll*min(Da,Db);
 			contactPhysics->momentRotationLaw=(sdec1->momentRotationLaw && sdec2->momentRotationLaw);
 			//contactPhysics->elasticRollingLimit = elasticRollingLimit;
 		}

=== modified file 'pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.hpp'
--- pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.hpp	2010-12-03 12:00:34 +0000
+++ pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.hpp	2010-12-08 10:43:28 +0000
@@ -21,11 +21,7 @@
 		YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_2xCohFrictMat_CohFrictPhys,IPhysFunctor,
 		"Generates cohesive-frictional interactions with moments. Used in the contact law :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment`.",
 		((bool,setCohesionNow,false,,"If true, assign cohesion to all existing contacts in current time-step. The flag is turned false automatically, so that assignment is done in the current timestep only."))
-		((bool,setCohesionOnNewContacts,false,,"If true, assign cohesion at all new contacts. If false, only existing contacts can be cohesive (also see :yref:`Ip2_2xCohFrictMat_CohFrictPhys::setCohesionNow`), and new contacts are only frictional."))
-		((Real,alphaKr,2.0,,"Dimensionless coefficient used for the rolling stiffness."))
-		((Real,alphaKtw,2.0,,"Dimensionless coefficient used for the twist stiffness."))
-		((Real,etaRoll,-1.,,"Dimensionless coefficient used to calculate the plastic rolling moment (if negative, plasticity will not be applied)."))
-		((bool,useMeanRad,false,,"use $meanRad^2$ to define rotational stifnesses, else use R1*R2"))
+		((bool,setCohesionOnNewContacts,false,,"If true, assign cohesion at all new contacts. If false, only existing contacts can be cohesive (also see :yref:`Ip2_2xCohFrictMat_CohFrictPhys::setCohesionNow`), and new contacts are only frictional."))	
 		,
 		cohesionDefinitionIteration = -1;
 		);