yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01320
[svn] r1807 - trunk/pkg/dem/Engine/EngineUnit
Author: richefeu
Date: 2009-06-23 16:13:28 +0200 (Tue, 23 Jun 2009)
New Revision: 1807
Modified:
trunk/pkg/dem/Engine/EngineUnit/BasicViscoelasticRelationships.cpp
Log:
Use the clump mass instead of clump-member mass for the determination of the effective mass.
Modified: trunk/pkg/dem/Engine/EngineUnit/BasicViscoelasticRelationships.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/BasicViscoelasticRelationships.cpp 2009-06-23 07:40:04 UTC (rev 1806)
+++ trunk/pkg/dem/Engine/EngineUnit/BasicViscoelasticRelationships.cpp 2009-06-23 14:13:28 UTC (rev 1807)
@@ -37,16 +37,45 @@
ViscoelasticInteraction* contactPhysics = YADE_CAST<ViscoelasticInteraction*>(interaction->interactionPhysics.get());
// Check that cn is dimensionless
- if (sdec1->cn >= 1.0) std::cout << "Warning: non dimensionless value for cn" << std::endl;
+ if (sdec1->cn >= 1.0 || sdec2->cn >= 1.0) std::cout << "Warning: non dimensionless value for cn" << std::endl;
- // Arbitrare ponderation...
+ // FIXME Arbitrary ponderations... (Should be replaced by parameters that depends on groupMasks)
contactPhysics->kn = 0.5 * (sdec1->kn + sdec2->kn);
contactPhysics->ks = 0.5 * (sdec1->ks + sdec2->ks);
contactPhysics->cn = 0.5 * (sdec1->cn + sdec2->cn);
contactPhysics->cs = 0.0; // not viscosity in the tangent direction
- // in fact, dimensionless cn is a ponderation of the critical value cn_crit = 2sqrt(kn*meff)
- contactPhysics->cn *= 2.0 * sqrt(contactPhysics->kn * ((sdec1->mass * sdec2->mass) / (sdec1->mass + sdec2->mass)));
+ // in fact, at this stage cn is a ponderation of the critical value cn_crit = 2sqrt(kn*meff)
+
+ shared_ptr<Body> bdy1 = (*((Omega::instance().getRootBody().get())->bodies))[ interaction->getId1() ];
+ shared_ptr<Body> bdy2 = (*((Omega::instance().getRootBody().get())->bodies))[ interaction->getId2() ];
+ double m1,m2;
+
+ if (bdy1->isClumpMember())
+ {
+ const shared_ptr<Body>& clump = (*((Omega::instance().getRootBody().get())->bodies))[ bdy1->clumpId ];
+ RigidBodyParameters* clumpRBP=YADE_CAST<RigidBodyParameters*> ( clump->physicalParameters.get() );
+ m1 = clumpRBP->mass;
+ }
+ else
+ { m1 = sdec1->mass ; }
+
+ if (bdy2->isClumpMember())
+ {
+ const shared_ptr<Body>& clump = (*((Omega::instance().getRootBody().get())->bodies))[ bdy2->clumpId ];
+ RigidBodyParameters* clumpRBP=YADE_CAST<RigidBodyParameters*> ( clump->physicalParameters.get() );
+ m2 = clumpRBP->mass;
+ }
+ else
+ { m2 = sdec2->mass; }
+
+ if (!bdy1->isDynamic) m1 *= 1e6;
+ if (!bdy2->isDynamic) m2 *= 1e6;
+
+ //cout << "m1 = " << m1 << endl;
+ //cout << "m2 = " << m2 << endl;
+
+ contactPhysics->cn *= 2.0 * sqrt(contactPhysics->kn * ((m1 * m2) / (m1 + m2)));
contactPhysics->tangensOfFrictionAngle = std::tan(std::min(sdec1->frictionAngle, sdec2->frictionAngle));
