yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2077
committer: Luc Scholtes <sch50p@fluent-ph>
branch nick: trunk
timestamp: Fri 2010-03-12 17:16:16 +1000
message:
  correction of few bugs in CohesiveFrictionalPM
modified:
  pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp


--
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/Engine/GlobalEngine/CohesiveFrictionalPM.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp	2010-03-09 23:43:53 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp	2010-03-12 07:16:16 +0000
@@ -22,14 +22,13 @@
 	Real dt = Omega::instance().getTimeStep();
 	
 	Real displN = geom->penetrationDepth; // NOTE: the sign for penetrationdepth is different from ScGeom and Dem3DofGeom: geom->penetrationDepth>0 when spheres interpenetrate
-	Real D = 0; // initial interparticular distance
 	Real Dtensile=phys->FnMax/phys->kn;
 	Real Dsoftening = phys->strengthSoftening*Dtensile; 
 	
 	/*to set the equilibrium distance between all cohesive elements when they first meet*/
-	if ( contact->isFresh(rootBody) && phys->isCohesive) { phys->initD = displN; phys->normalForce = Vector3r::ZERO; phys->shearForce = Vector3r::ZERO;}
+	if ( contact->isFresh(rootBody) ) { phys->initD = displN; phys->normalForce = Vector3r::ZERO; phys->shearForce = Vector3r::ZERO;}
 	
-	D = displN - phys->initD; // displacement is computed depending on the equilibrium distance
+	Real D = displN - phys->initD; // interparticular distance is computed depending on the equilibrium distance (which is 0 for non cohesive interactions)
 	
 	/* Determination of interaction */
 	if (D < 0){ //spheres do not touch 
@@ -45,7 +44,7 @@
 		
 	if ((D < 0) && (abs(D) > Dtensile)) { //to take into account strength softening
 	  Dsoft = D+Dtensile; // Dsoft<0 for a negative value of Fn (attractive force)
-	  Fn = phys->FnMax+(phys->kn/phys->strengthSoftening)*Dsoft; // computes FnMax - FnSoftening
+	  Fn = -(phys->FnMax+(phys->kn/phys->strengthSoftening)*Dsoft); // computes FnMax - FnSoftening
 	}
 	else {
 	  Fn = phys->kn*D;
@@ -69,8 +68,7 @@
 	Real maxFs = abs(Fn)*phys->tanFrictionAngle + phys->FsMax;
 		
 	if (Fs.SquaredLength() > maxFs*maxFs){ 
-	  phys->isCohesive=false; //if not already false to delete cohesive interactions
-	  Fs*=maxFs/Fs.Length(); // to fix the shear force to its yielding value with the right sign
+	  Fs*=maxFs/Fs.Length(); // to fix the shear force to its yielding value
 	}
 	phys->shearForce = Fs;