yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2105
committer: Sergei D. <sega@think>
branch nick: trunk
timestamp: Thu 2010-03-25 23:47:48 +0300
message:
  yet more fixing ViscoelasticPM
modified:
  pkg/dem/meta/ViscoelasticPM.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/meta/ViscoelasticPM.cpp'
--- pkg/dem/meta/ViscoelasticPM.cpp	2010-03-25 14:55:30 +0000
+++ pkg/dem/meta/ViscoelasticPM.cpp	2010-03-25 20:47:48 +0000
@@ -77,21 +77,28 @@
 	// As Chiara Modenese suggest, we store the elastic part 
 	// and then add the viscous part if we pass the Mohr-Coulomb criterion.
 	// See http://www.mail-archive.com/yade-users@xxxxxxxxxxxxxxxxxxx/msg01391.html
-	shearForce += phys.ks*dt*shearVelocity;
-	Vector3r shearForceVisc = Vector3r::ZERO;
+	shearForce += phys.ks*dt*shearVelocity; // the elastic shear force have a history, but
+	Vector3r shearForceVisc = Vector3r::ZERO; // the viscous shear damping haven't a history because it is a function of the instant velocity 
 
 	phys.normalForce = ( phys.kn * geom.penetrationDepth + phys.cn * normalVelocity ) * geom.normal;
 	phys.prevNormal = geom.normal;
 
-	Real maxFs = phys.normalForce.SquaredLength() * std::pow(phys.tangensOfFrictionAngle,2);
+	const Real maxFs = phys.normalForce.SquaredLength() * std::pow(phys.tangensOfFrictionAngle,2);
 	if( shearForce.SquaredLength() > maxFs )
 	{
-		maxFs = Mathr::Sqrt(maxFs) / shearForce.Length();
-		shearForce *= maxFs;
-		shearForceVisc =  phys.cs*shearVelocity;
+		// Then Mohr-Coulomb is violated (so, we slip), 
+		// we have the max value of the shear force, so 
+		// we consider only friction damping.
+		const Real ratio = Mathr::Sqrt(maxFs) / shearForce.Length();
+		shearForce *= ratio;
+	} 
+	else 
+	{
+		// Then no slip occurs we consider friction damping + viscous damping.
+		shearForceVisc = phys.cs*shearVelocity; 
 	}
 
-	Vector3r f = phys.normalForce + shearForce + shearForceVisc;
+	const Vector3r f = phys.normalForce + shearForce + shearForceVisc;
 	addForce (id1,-f,rootBody);
 	addForce (id2, f,rootBody);
 	addTorque(id1,-c1x.Cross(f),rootBody);