Author: eudoxos
Date: 2009-03-09 21:57:46 +0100 (Mon, 09 Mar 2009)
New Revision: 1715
Modified:
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
trunk/pkg/dem/SConscript
Log:
1. Added SpheresContactGeometry::updateShearForce, will be used (not activated
though yet) by ElasticContactLaw and other.
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
2009-03-08 23:03:42 UTC (rev 1714)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
2009-03-09 20:57:46 UTC (rev 1715)
@@ -9,6 +9,62 @@
// At least one virtual method must be in the .cpp file (!!!)
SpheresContactGeometry::~SpheresContactGeometry(){};
+void SpheresContactGeometry::updateShearForce(Vector3r& shearForce, Real ks,
const Vector3r& prevNormal, const RigidBodyParameters* rbp1, const
RigidBodyParameters* rbp2, bool isDynamic1, bool isDynamic2, Real dt, bool
avoidGranularRatcheting){
+
+ Vector3r axis;
+ Real angle;
+
+ // approximated rotations
+ axis = prevNormal.Cross(normal);
+ shearForce -= shearForce.Cross(axis);
+ //angle =
dt*0.5*currentContactGeometry->normal.Dot(de1->angularVelocity+de2->angularVelocity);
+ //FIXME: if one body is kinematic then assumed its rotation
centre does not lies along the normal
+ //(i.e. virtual sphere, which replaces this kinematic body on
contact, does not rotate)
+ Vector3r summaryAngularVelocity(Vector3r::ZERO);
+ if (isDynamic1) summaryAngularVelocity += rbp1->angularVelocity;
+ if (isDynamic2) summaryAngularVelocity += rbp2->angularVelocity;
+ angle = dt*0.5*normal.Dot(summaryAngularVelocity);
+ axis = angle*normal;
+ shearForce -= shearForce.Cross(axis);
+
+ // exact rotations
+ #if 0
+ Quaternionr q;
+ axis =
prevNormal.Cross(normal);
+ angle =
acos(normal.Dot(prevNormal));
+ q.FromAngleAxis(angle,axis);
+ shearForce = shearForce*q;
+ angle =
dt*0.5*normal.dot(rbp1->angularVelocity+rbp2->angularVelocity);
+ axis = normal;
+ q.FromAngleAxis(angle,axis);
+ shearForce = q*shearForce;
+ #endif
+
+ Vector3r& x = contactPoint;
+ Vector3r c1x, c2x;
+
+ if(avoidGranularRatcheting){
+ /* The following definition of c1x and c2x is to avoid
"granular ratcheting"
+ * (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN,
+ * "Micro-mechanical investigation of granular ratcheting, in
Cyclic Behaviour of Soils and Liquefaction Phenomena",
+ * ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 -
and a lot more papers from the same authors) */
+ c1x = (isDynamic1) ? radius1*normal : x - rbp1->zeroPoint;
+ c2x = (isDynamic2) ? -radius2*normal : x - rbp2->zeroPoint;
+ }
+ else {
+ c1x = (x - rbp1->se3.position);
+ c2x = (x - rbp2->se3.position);
+ }
+
+ Vector3r relativeVelocity =
(rbp2->velocity+rbp2->angularVelocity.Cross(c2x)) -
(rbp1->velocity+rbp1->angularVelocity.Cross(c1x));
+ Vector3r shearVelocity =
relativeVelocity-normal.Dot(relativeVelocity)*normal;
+ Vector3r shearDisplacement = shearVelocity*dt;
+ shearForce -= ks*shearDisplacement;
+}
+
+
+
+
Vector3r SpheresContactGeometry::relRotVector() const{
Quaternionr relOri12=ori1.Conjugate()*ori2;
Quaternionr oriDiff=initRelOri12.Conjugate()*relOri12;
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
2009-03-08 23:03:42 UTC (rev 1714)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
2009-03-09 20:57:46 UTC (rev 1715)
@@ -5,6 +5,7 @@
#include<yade/core/InteractionGeometry.hpp>
#include<yade/lib-base/yadeWm3.hpp>
+#include<yade/pkg-common/RigidBodyParameters.hpp>
/*! Class representing geometry of two spheres in contact.
*
* exactRot code
@@ -98,6 +99,9 @@
SpheresContactGeometry():contactPoint(Vector3r::ZERO),radius1(0),radius2(0),facetContactFace(0.),hasShear(false),pos1(Vector3r::ZERO),pos2(Vector3r::ZERO),ori1(Quaternionr::IDENTITY),ori2(Quaternionr::IDENTITY),cp1rel(Quaternionr::IDENTITY),cp2rel(Quaternionr::IDENTITY),d1(0),d2(0),d0(0),initRelOri12(Quaternionr::IDENTITY){createIndex();}
virtual ~SpheresContactGeometry();
+
+ void updateShearForce(Vector3r& shearForce, Real ks, const
Vector3r& prevNormal, const RigidBodyParameters* rbp1, const
RigidBodyParameters* rbp2, bool isDynamic1, bool isDynamic2, Real dt, bool
avoidGranularRatcheting=true);
+
REGISTER_ATTRIBUTES(/* no attributes from base class */,
(normal)
(contactPoint)
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-03-08
23:03:42 UTC (rev 1714)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-03-09
20:57:46 UTC (rev 1715)
@@ -114,6 +114,10 @@
Real un=currentContactGeometry->penetrationDepth;
currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real)
0)*currentContactGeometry->normal;
+ #if 0
+ // the core under #else, refactored
+
currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,isDynamic1,isDynamic2,dt);
+ #else
Vector3r axis;
Real angle;
@@ -163,6 +167,8 @@
Vector3r shearVelocity =
relativeVelocity-currentContactGeometry->normal.Dot(relativeVelocity)*currentContactGeometry->normal;
Vector3r shearDisplacement =
shearVelocity*dt;
shearForce -=
currentContactPhysics->ks*shearDisplacement;
+
+ #endif
// PFC3d SlipModel, is using friction angle. CoulombCriterion
Real maxFs =
currentContactPhysics->normalForce.SquaredLength() *
std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
@@ -173,7 +179,9 @@
}
////////// PFC3d SlipModel
- Vector3r f =
currentContactPhysics->normalForce + shearForce;
+ Vector3r f=currentContactPhysics->normalForce +
shearForce;
+ Vector3r
c1x(currentContactGeometry->contactPoint-de1->se3.position),
+
c2x(currentContactGeometry->contactPoint-de2->se3.position);
#ifdef BEX_CONTAINER
ncb->bex.addForce (id1,-f);
ncb->bex.addForce (id2,+f);
Modified: trunk/pkg/dem/SConscript
===================================================================
--- trunk/pkg/dem/SConscript 2009-03-08 23:03:42 UTC (rev 1714)
+++ trunk/pkg/dem/SConscript 2009-03-09 20:57:46 UTC (rev 1715)
@@ -54,15 +54,15 @@
env.SharedLibrary('SpheresContactGeometry',
['DataClass/InteractionGeometry/SpheresContactGeometry.cpp'],
- LIBS=env['LIBS']+['yade-serialization','yade-base']),
+ LIBS=env['LIBS']+['RigidBodyParameters']),
env.SharedLibrary('SDECLinkGeometry',
- ['DataClass/InteractionGeometry/SDECLinkGeometry.cpp'],
- LIBS=env['LIBS']+['yade-serialization','yade-base']),
+ ['DataClass/InteractionGeometry/SDECLinkGeometry.cpp']
+ ),
env.SharedLibrary('InteractionOfMyTetrahedron',
['DataClass/InteractionGeometry/InteractionOfMyTetrahedron.cpp'],
- LIBS=env['LIBS']+['yade-serialization', 'yade-base']),
+ ),
env.SharedLibrary('ElasticContactInteraction',
['DataClass/InteractionPhysics/ElasticContactInteraction.cpp'],
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
