Hi, all: Attached are two files for viscous damping, I will make some documentation later and I will make it compatible with latest SVN when have the svn access.
Feng Chen Graduate Student Department of Civil and Environmental Engineering 223 Perkins Hall University of Tennessee, Knoxville, 37996 http://fchen3.googlepages.com/home
/************************************************************************* * Copyright (C) 2008 by Feng Chen ([email protected]) * * Department of Civil and Environmental Engineering * * 223 Perkins Hall * * University of Tennessee, Knoxville, 37996 * * http://fchen3.googlepages.com/home * * * * This program is free software; it is licensed under the terms of the * * GNU General Public License v2 or later. See file LICENSE for details. * *************************************************************************/ #ifndef VISCOUS_FORCE_DAMPING_HPP #define VISCOUS_FORCE_DAMPING_HPP #include<yade/core/InteractionSolver.hpp> #include <set> #include <boost/tuple/tuple.hpp> class PhysicalAction; class ViscousForceDamping : public InteractionSolver { /// Attributes private : shared_ptr<PhysicalAction> actionForce; shared_ptr<PhysicalAction> actionMomentum; public : Real betaNormal; Real betaShear; int sdecGroupMask; bool momentRotationLaw; ViscousForceDamping(); void action(Body* body); protected : void registerAttributes(); REGISTER_CLASS_NAME(ViscousForceDamping); REGISTER_BASE_CLASS_NAME(InteractionSolver); }; REGISTER_SERIALIZABLE(ViscousForceDamping,false); #endif // VISCOUS_FORCE_DAMPING_HPP
/************************************************************************* * Copyright (C) 2008 by Feng Chen ([email protected]) * * Department of Civil and Environmental Engineering * * 223 Perkins Hall * * University of Tennessee, Knoxville, 37996 * * http://fchen3.googlepages.com/home * * * * This program is free software; it is licensed under the terms of the * * GNU General Public License v2 or later. See file LICENSE for details. * *************************************************************************/ #include "ViscousForceDamping.hpp" #include<yade/pkg-dem/BodyMacroParameters.hpp> #include<yade/pkg-dem/SpheresContactGeometry.hpp> #include<yade/pkg-dem/SDECLinkGeometry.hpp> #include<yade/pkg-dem/ElasticContactInteraction.hpp> #include<yade/pkg-dem/SDECLinkPhysics.hpp> #include<yade/core/Omega.hpp> #include<yade/core/MetaBody.hpp> #include<yade/pkg-common/Force.hpp> #include<yade/pkg-common/Momentum.hpp> #include<yade/core/PhysicalAction.hpp> ViscousForceDamping::ViscousForceDamping() : InteractionSolver() , actionForce(new Force) , actionMomentum(new Momentum), betaNormal(0.0), betaShear(0.0) { sdecGroupMask=1; momentRotationLaw = true; ///==================================== ///betaNormal = 0.0; ///betaShear = 0.0; } void ViscousForceDamping::registerAttributes() { InteractionSolver::registerAttributes(); REGISTER_ATTRIBUTE(sdecGroupMask); REGISTER_ATTRIBUTE(momentRotationLaw); REGISTER_ATTRIBUTE(betaNormal); REGISTER_ATTRIBUTE(betaShear); } //FIXME : remove bool first !!!!! void ViscousForceDamping::action(Body* body) { MetaBody * ncb = YADE_CAST<MetaBody*>(body); shared_ptr<BodyContainer>& bodies = ncb->bodies; Real dt = Omega::instance().getTimeStep(); /// Non Permanents Links /// InteractionContainer::iterator ii = ncb->transientInteractions->begin(); InteractionContainer::iterator iiEnd = ncb->transientInteractions->end(); for( ; ii!=iiEnd ; ++ii ) { if ((*ii)->isReal) { const shared_ptr<Interaction>& contact = *ii; int id1 = contact->getId1(); int id2 = contact->getId2(); if( !( (*bodies)[id1]->getGroupMask() & (*bodies)[id2]->getGroupMask() & sdecGroupMask) ) continue; // skip other groups, BTW: this is example of a good usage of 'continue' keyword BodyMacroParameters* de1 = YADE_CAST<BodyMacroParameters*>((*bodies)[id1]->physicalParameters.get()); BodyMacroParameters* de2 = YADE_CAST<BodyMacroParameters*>((*bodies)[id2]->physicalParameters.get()); SpheresContactGeometry* currentContactGeometry = YADE_CAST<SpheresContactGeometry*>(contact->interactionGeometry.get()); ElasticContactInteraction* currentContactPhysics = YADE_CAST<ElasticContactInteraction*> (contact->interactionPhysics.get()); Vector3r& shearForce = currentContactPhysics->shearForce; if ( contact->isNew) shearForce = Vector3r(0,0,0); Real un = currentContactGeometry->penetrationDepth; currentContactPhysics->normalForce = currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal; Vector3r axis; Real angle; /// Here is the code with approximated rotations /// axis = currentContactPhysics->prevNormal.Cross(currentContactGeometry->normal); shearForce -= shearForce.Cross(axis); angle = dt*0.5*currentContactGeometry->normal.Dot(de1->angularVelocity+de2->angularVelocity); axis = angle*currentContactGeometry->normal; shearForce -= shearForce.Cross(axis); /// Here is the code with exact rotations /// // Quaternionr q; // // axis = currentContactPhysics->prevNormal.cross(currentContactGeometry->normal); // angle = acos(currentContactGeometry->normal.dot(currentContactPhysics->prevNormal)); // q.fromAngleAxis(angle,axis); // // currentContactPhysics->shearForce = currentContactPhysics->shearForce*q; // // angle = dt*0.5*currentContactGeometry->normal.dot(de1->angularVelocity+de2->angularVelocity); // axis = currentContactGeometry->normal; // q.fromAngleAxis(angle,axis); // currentContactPhysics->shearForce = q*currentContactPhysics->shearForce; /// /// Vector3r x = currentContactGeometry->contactPoint; Vector3r c1x = (x - de1->se3.position); Vector3r c2x = (x - de2->se3.position); Vector3r relativeVelocity = (de2->velocity+de2->angularVelocity.Cross(c2x)) - (de1->velocity+de1->angularVelocity.Cross(c1x)); Vector3r shearVelocity = relativeVelocity-currentContactGeometry->normal.Dot(relativeVelocity)*currentContactGeometry->normal; Vector3r shearDisplacement = shearVelocity*dt; shearForce -= currentContactPhysics->ks*shearDisplacement; //=============================================================================== //Viscous damping Vector3r normalVelocity = relativeVelocity - shearVelocity; Real m1 = de1->mass; Real m2 = de2->mass; Real mBar = m1*m2/(m1+m2); if (!(*bodies)[id1]->isDynamic) ///id1 is a static wall mBar = m2; if (!(*bodies)[id2]->isDynamic) ///id2 is a static wall mBar = m1; // std::cerr << "kn:" << currentContactPhysics->kn << endl; // std::cerr << "ks:" << currentContactPhysics->ks << endl; // std::cerr << "mBar:" << mBar << endl; Real cCritNormal = 2.0*sqrt(mBar*currentContactPhysics->kn); Real cCritShear = 2.0*sqrt(mBar*currentContactPhysics->ks); Real cNormal = betaNormal*cCritNormal; Real cShear = betaShear*cCritShear; Vector3r normalDampingForce = -cNormal*normalVelocity; Vector3r shearDampingForce = -cShear*shearVelocity; Vector3r viscousDampingForce = normalDampingForce + shearDampingForce; // Add forces static_cast<Force*> ( ncb->physicalActions->find( id1 , actionForce ->getClassIndex() ).get() )->force -= viscousDampingForce; static_cast<Force*> ( ncb->physicalActions->find( id2 , actionForce ->getClassIndex() ).get() )->force += viscousDampingForce; static_cast<Momentum*>( ncb->physicalActions->find( id1 , actionMomentum->getClassIndex() ).get() )->momentum -= c1x.Cross(viscousDampingForce); static_cast<Momentum*>( ncb->physicalActions->find( id2 , actionMomentum->getClassIndex() ).get() )->momentum += c2x.Cross(viscousDampingForce); currentContactPhysics->prevNormal = currentContactGeometry->normal; } } }
_______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
