------------------------------------------------------------ revno: 3804 committer: Klaus Thoeni <[email protected]> timestamp: Sat 2014-01-25 08:59:41 +1100 message: new simple contact law with normal viscose damping which allows to specify kn and ks/kn in Ip2 added: pkg/dem/FrictViscoPM.cpp pkg/dem/FrictViscoPM.hpp
-- lp:yade https://code.launchpad.net/~yade-pkg/yade/git-trunk Your team Yade developers is subscribed to branch lp:yade. To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== added file 'pkg/dem/FrictViscoPM.cpp' --- pkg/dem/FrictViscoPM.cpp 1970-01-01 00:00:00 +0000 +++ pkg/dem/FrictViscoPM.cpp 2014-01-24 21:59:41 +0000 @@ -0,0 +1,228 @@ +/************************************************************************* +* Copyright (C) 2014 by Klaus Thoeni * +* [email protected] * +* * +* 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"FrictViscoPM.hpp" +#include<yade/pkg/dem/ScGeom.hpp> +#include<yade/core/Omega.hpp> +#include<yade/core/Scene.hpp> + +YADE_PLUGIN((FrictViscoMat)(FrictViscoPhys)(Ip2_FrictViscoMat_FrictViscoMat_FrictViscoPhys)(Ip2_FrictMat_FrictViscoMat_FrictViscoPhys)(Law2_ScGeom_FrictViscoPhys_CundallStrackVisco)); + +FrictViscoMat::~FrictViscoMat(){} + +/********************** Ip2_FrictViscoMat_FrictMat_FrictViscoPhys ****************************/ +CREATE_LOGGER(FrictViscoPhys); + +FrictViscoPhys::~FrictViscoPhys(){}; + +/********************** Ip2_FrictViscoMat_FrictMat_FrictViscoPhys ****************************/ +CREATE_LOGGER(Ip2_FrictViscoMat_FrictViscoMat_FrictViscoPhys); + +void Ip2_FrictViscoMat_FrictViscoMat_FrictViscoPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction){ + + LOG_TRACE( "Ip2_FrictViscoMat_FrictViscoMat_FrictViscoPhys::go - contact law" ); + + if(interaction->phys) return; + const shared_ptr<FrictViscoMat>& mat1 = YADE_PTR_CAST<FrictViscoMat>(b1); + const shared_ptr<FrictViscoMat>& mat2 = YADE_PTR_CAST<FrictViscoMat>(b2); + interaction->phys = shared_ptr<FrictViscoPhys>(new FrictViscoPhys()); + const shared_ptr<FrictViscoPhys>& contactPhysics = YADE_PTR_CAST<FrictViscoPhys>(interaction->phys); + Real Ea = mat1->young; + Real Eb = mat2->young; + Real Va = mat1->poisson; + Real Vb = mat2->poisson; + + Real Ra,Rb; + assert(dynamic_cast<GenericSpheresContact*>(interaction->geom.get()));//only in debug mode + GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->geom.get()); + Ra=sphCont->refR1>0?sphCont->refR1:sphCont->refR2; + Rb=sphCont->refR2>0?sphCont->refR2:sphCont->refR1; + + // calculate stiffness from MatchMaker or use harmonic avarage as usual if not given + Real Kn = (kn) ? (*kn)(mat1->id,mat2->id) : 2.*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb); + Real Ks = (kRatio) ? (*kRatio)(mat1->id,mat2->id)*Kn : 2.*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb); + + Real frictionAngle = (!frictAngle) ? std::min(mat1->frictionAngle,mat2->frictionAngle) : (*frictAngle)(mat1->id,mat2->id,mat1->frictionAngle,mat2->frictionAngle); + contactPhysics->tangensOfFrictionAngle = std::tan(frictionAngle); + contactPhysics->kn = Kn; + contactPhysics->ks = Ks; + + /************************/ + /* DAMPING COEFFICIENTS */ + /************************/ + Real betana = mat1->betan; + Real betanb = mat2->betan; + + // inclusion of local viscous damping + if (betana!=0 || betanb!=0){ + Body::id_t ida = interaction->getId1(); // get id body 1 + Body::id_t idb = interaction->getId2(); // get id body 2 + State* dea = Body::byId(ida,scene)->state.get(); + State* deb = Body::byId(idb,scene)->state.get(); + const shared_ptr<Body>& ba=Body::byId(ida,scene); + const shared_ptr<Body>& bb=Body::byId(idb,scene); + Real mbar = (!ba->isDynamic() && bb->isDynamic()) ? deb->mass : ((!bb->isDynamic() && ba->isDynamic()) ? dea->mass : (dea->mass*deb->mass / (dea->mass + deb->mass))); // get equivalent mass if both bodies are dynamic, if not set it equal to the one of the dynamic body + TRVAR2(Kn,mbar); + contactPhysics->cn_crit = 2.*sqrt(mbar*Kn); // Critical damping coefficient (normal direction) + contactPhysics->cn = contactPhysics->cn_crit * ( (betana!=0 && betanb!=0) ? ((betana+betanb)/2.) : ( (betanb==0) ? betana : betanb )); // Damping normal coefficient + } + else + contactPhysics->cn=0.; + TRVAR1(contactPhysics->cn); + +} + +/********************** Ip2_FrictViscoMat_FrictMat_FrictViscoPhys ****************************/ +CREATE_LOGGER(Ip2_FrictMat_FrictViscoMat_FrictViscoPhys); + +void Ip2_FrictMat_FrictViscoMat_FrictViscoPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction){ + + LOG_TRACE( "Ip2_FrictMat_FrictViscoMat_FrictViscoPhys::go - contact law" ); + + if(interaction->phys) return; + const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1); + const shared_ptr<FrictViscoMat>& mat2 = YADE_PTR_CAST<FrictViscoMat>(b2); + interaction->phys = shared_ptr<FrictViscoPhys>(new FrictViscoPhys()); + const shared_ptr<FrictViscoPhys>& contactPhysics = YADE_PTR_CAST<FrictViscoPhys>(interaction->phys); + Real Ea = mat1->young; + Real Eb = mat2->young; + Real Va = mat1->poisson; + Real Vb = mat2->poisson; + + Real Ra,Rb; + assert(dynamic_cast<GenericSpheresContact*>(interaction->geom.get()));//only in debug mode + GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->geom.get()); + Ra=sphCont->refR1>0?sphCont->refR1:sphCont->refR2; + Rb=sphCont->refR2>0?sphCont->refR2:sphCont->refR1; + + // calculate stiffness from MatchMaker or use harmonic avarage as usual if not given + Real Kn = (kn) ? (*kn)(mat1->id,mat2->id) : 2.*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb); + Real Ks = (kRatio) ? (*kRatio)(mat1->id,mat2->id)*Kn : 2.*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb); + + Real frictionAngle = (!frictAngle) ? std::min(mat1->frictionAngle,mat2->frictionAngle) : (*frictAngle)(mat1->id,mat2->id,mat1->frictionAngle,mat2->frictionAngle); + contactPhysics->tangensOfFrictionAngle = std::tan(frictionAngle); + contactPhysics->kn = Kn; + contactPhysics->ks = Ks; + + /************************/ + /* DAMPING COEFFICIENTS */ + /************************/ + Real betanb = mat2->betan; + + // inclusion of local viscous damping + if (betanb!=0){ + Body::id_t ida = interaction->getId1(); // get id body 1 + Body::id_t idb = interaction->getId2(); // get id body 2 + State* dea = Body::byId(ida,scene)->state.get(); + State* deb = Body::byId(idb,scene)->state.get(); + const shared_ptr<Body>& ba=Body::byId(ida,scene); + const shared_ptr<Body>& bb=Body::byId(idb,scene); + Real mbar = (!ba->isDynamic() && bb->isDynamic()) ? deb->mass : ((!bb->isDynamic() && ba->isDynamic()) ? dea->mass : (dea->mass*deb->mass / (dea->mass + deb->mass))); // get equivalent mass if both bodies are dynamic, if not set it equal to the one of the dynamic body + TRVAR2(Kn,mbar); + contactPhysics->cn_crit = 2.*sqrt(mbar*Kn); // Critical damping coefficient (normal direction) + contactPhysics->cn = contactPhysics->cn_crit * betanb; // Damping normal coefficient + } + else + contactPhysics->cn=0.; + TRVAR1(contactPhysics->cn); + +} + +/********************** Law2_ScGeom_FrictViscoPhys_CundallStrackVisco ****************************/ + +// #if 1 +Real Law2_ScGeom_FrictViscoPhys_CundallStrackVisco::getPlasticDissipation() {return (Real) plasticDissipation;} +void Law2_ScGeom_FrictViscoPhys_CundallStrackVisco::initPlasticDissipation(Real initVal) {plasticDissipation.reset(); plasticDissipation+=initVal;} +Real Law2_ScGeom_FrictViscoPhys_CundallStrackVisco::elasticEnergy() +{ + Real energy=0; + FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){ + if(!I->isReal()) continue; + FrictPhys* phys = dynamic_cast<FrictPhys*>(I->phys.get()); + if(phys) { + energy += 0.5*(phys->normalForce.squaredNorm()/phys->kn + phys->shearForce.squaredNorm()/phys->ks);} + } + return energy; +} +// #endif + + +CREATE_LOGGER(Law2_ScGeom_FrictViscoPhys_CundallStrackVisco); + +void Law2_ScGeom_FrictViscoPhys_CundallStrackVisco::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact) +{ + LOG_TRACE( "Law2_ScGeom_FrictViscoPhys_CundallStrackVisco::go - create interaction physics" ); + + int id1 = contact->getId1(), id2 = contact->getId2(); + + ScGeom* geom= static_cast<ScGeom*>(ig.get()); + FrictViscoPhys* phys = static_cast<FrictViscoPhys*>(ip.get()); + if(geom->penetrationDepth <0){ + if (neverErase) { + phys->shearForce = Vector3r::Zero(); + phys->normalForce = Vector3r::Zero();} + else scene->interactions->requestErase(contact); + return; + } + Real& un=geom->penetrationDepth; + phys->normalForce = phys->kn*std::max(un,(Real) 0) * geom->normal; + + /************************/ + /* DAMPING CONTRIBUTION */ + /************************/ + + // define shifts to handle periodicity + const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero(); + const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(contact->cellDist): Vector3r::Zero(); + + State* de1 = Body::byId(id1,scene)->state.get(); + State* de2 = Body::byId(id2,scene)->state.get(); + + // get incident velocity + Vector3r incidentV = geom->getIncidentVel(de1, de2, scene->dt, shift2, shiftVel); + Vector3r incidentVn = geom->normal.dot(incidentV)*geom->normal; // contact normal velocity + phys->normalViscous = phys->cn*incidentVn; + phys->normalForce -= phys->normalViscous; + + // shear force + Vector3r& shearForce = geom->rotate(phys->shearForce); + const Vector3r& shearDisp = geom->shearIncrement(); + shearForce -= phys->ks*shearDisp; + Real maxFs = phys->normalForce.squaredNorm()*std::pow(phys->tangensOfFrictionAngle,2); + + if (!scene->trackEnergy && !traceEnergy){//Update force but don't compute energy terms (see below)) + // PFC3d SlipModel, is using friction angle. CoulombCriterion + if( shearForce.squaredNorm() > maxFs ){ + Real ratio = sqrt(maxFs) / shearForce.norm(); + shearForce *= ratio;} + } else { + //almost the same with additional Vector3r instatinated for energy tracing, + //duplicated block to make sure there is no cost for the instanciation of the vector when traceEnergy==false + if(shearForce.squaredNorm() > maxFs){ + Real ratio = sqrt(maxFs) / shearForce.norm(); + Vector3r trialForce=shearForce;//store prev force for definition of plastic slip + //define the plastic work input and increment the total plastic energy dissipated + shearForce *= ratio; + Real dissip=((1/phys->ks)*(trialForce-shearForce))/*plastic disp*/ .dot(shearForce)/*active force*/; + if (traceEnergy) plasticDissipation += dissip; + else if(dissip>0) scene->energy->add(dissip,"plastDissip",plastDissipIx,/*reset*/false); + } + // compute elastic energy as well + scene->energy->add(0.5*(phys->normalForce.squaredNorm()/phys->kn+phys->shearForce.squaredNorm()/phys->ks),"elastPotential",elastPotentialIx,/*reset at every timestep*/true); + } + if (!scene->isPeriodic && !sphericalBodies) { + applyForceAtContactPoint(-phys->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);} + else {//we need to use correct branches in the periodic case, the following apply for spheres only + Vector3r force = -phys->normalForce-shearForce; + scene->forces.addForce(id1,force); + scene->forces.addForce(id2,-force); + scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force)); + scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force)); + } +} + === added file 'pkg/dem/FrictViscoPM.hpp' --- pkg/dem/FrictViscoPM.hpp 1970-01-01 00:00:00 +0000 +++ pkg/dem/FrictViscoPM.hpp 2014-01-24 21:59:41 +0000 @@ -0,0 +1,120 @@ +/************************************************************************* +* Copyright (C) 2014 by Klaus Thoeni * +* [email protected] * +* * +* 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. * +*************************************************************************/ + +/** +=== OVERVIEW OF FrictViscoPM === + +A particle model for friction and viscous damping in normal direction. The damping coefficient +has to be defined with the material whereas the contact stiffness kn and ks/kn can be defined with the Ip2 functor. + +Remarks: +- maybe there is a better way of implementing this without copying from ElasticContactLaw +- maybe we can combine some ideas of this contact law with other contact laws +*/ + +#pragma once + +#include<yade/pkg/common/ElastMat.hpp> +#include<yade/pkg/common/Dispatching.hpp> +#include<yade/pkg/common/MatchMaker.hpp> +#include<yade/pkg/dem/FrictPhys.hpp> +#include<yade/pkg/dem/ScGeom.hpp> +#include<yade/lib/base/openmp-accu.hpp> + +#include<set> +#include<boost/tuple/tuple.hpp> + +/** This class holds information associated with each body */ +class FrictViscoMat: public FrictMat { + public: + virtual ~FrictViscoMat(); + YADE_CLASS_BASE_DOC_ATTRS_CTOR(FrictViscoMat,FrictMat,"Material for use with the FrictViscoPM classes", + ((Real,betan,0.,,"Fraction of the viscous damping coefficient in normal direction equal to $\\frac{c_{n}}{C_{n,crit}}$.")) + , + createIndex(); + ); + DECLARE_LOGGER; + REGISTER_CLASS_INDEX(FrictViscoMat,FrictMat); +}; +REGISTER_SERIALIZABLE(FrictViscoMat); + + +/** This class holds information associated with each interaction */ +class FrictViscoPhys: public FrictPhys { + public: + virtual ~FrictViscoPhys(); + YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(FrictViscoPhys,FrictPhys,"Representation of a single interaction of the FrictViscoPM type, storage for relevant parameters", + ((Real,cn_crit,NaN,,"Normal viscous constant for ctitical damping defined as $\\c_{n}=C_{n,crit}\\beta_n$.")) + ((Real,cn,NaN,,"Normal viscous constant defined as $\\c_{n}=c_{n,crit}\\beta_n$.")) + ((Vector3r,normalViscous,Vector3r::Zero(),,"Normal viscous component")) + , + createIndex(); + , + ); + DECLARE_LOGGER; + REGISTER_CLASS_INDEX(FrictViscoPhys,FrictPhys); +}; +REGISTER_SERIALIZABLE(FrictViscoPhys); + +/** 2d functor creating IPhys (Ip2) taking FrictViscoMat and FrictViscoMat of 2 bodies, returning type FrictViscoPhys */ +class Ip2_FrictViscoMat_FrictViscoMat_FrictViscoPhys: public IPhysFunctor{ + public: + virtual void go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction); + + FUNCTOR2D(FrictViscoMat,FrictViscoMat); + + YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictViscoMat_FrictViscoMat_FrictViscoPhys,IPhysFunctor,"Converts 2 :yref:`FrictViscoMat` instances to :yref:`FrictViscoPhys` with corresponding parameters. Basically this functor corresponds to :yref:`Ip2_FrictMat_FrictMat_FrictPhys` with the only difference that damping in normal direction can be considered.", + ((shared_ptr<MatchMaker>,kn,,,"Instance of :yref:`MatchMaker` determining how to compute interaction's normal contact stiffnesses. If this value is not given the elastic properties (i.e. young) of the two colliding materials are used to calculate the stiffness.")) + ((shared_ptr<MatchMaker>,kRatio,,,"Instance of :yref:`MatchMaker` determining how to compute interaction's shear contact stiffnesses. If this value is not given the elastic properties (i.e. poisson) of the two colliding materials are used to calculate the stiffness.")) + ((shared_ptr<MatchMaker>,frictAngle,,,"Instance of :yref:`MatchMaker` determining how to compute interaction's friction angle. If ``None``, minimum value is used.")) + ); + DECLARE_LOGGER; +}; +REGISTER_SERIALIZABLE(Ip2_FrictViscoMat_FrictViscoMat_FrictViscoPhys); + +/** 2d functor creating IPhys (Ip2) taking FrictMat and FrictViscoMat of 2 bodies, returning type FrictViscoPhys */ +class Ip2_FrictMat_FrictViscoMat_FrictViscoPhys: public IPhysFunctor{ + public: + virtual void go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction); + + FUNCTOR2D(FrictMat,FrictViscoMat); + + YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictViscoMat_FrictViscoPhys,IPhysFunctor,"Converts a :yref:`FrictMat` and :yref:`FrictViscoMat` instance to :yref:`FrictViscoPhys` with corresponding parameters. Basically this functor corresponds to :yref:`Ip2_FrictMat_FrictMat_FrictPhys` with the only difference that damping in normal direction can be considered.", + ((shared_ptr<MatchMaker>,kn,,,"Instance of :yref:`MatchMaker` determining how to compute interaction's normal contact stiffnesses. If this value is not given the elastic properties (i.e. young) of the two colliding materials are used to calculate the stiffness.")) + ((shared_ptr<MatchMaker>,kRatio,,,"Instance of :yref:`MatchMaker` determining how to compute interaction's shear contact stiffnesses. If this value is not given the elastic properties (i.e. poisson) of the two colliding materials are used to calculate the stiffness.")) + ((shared_ptr<MatchMaker>,frictAngle,,,"Instance of :yref:`MatchMaker` determining how to compute interaction's friction angle. If ``None``, minimum value is used.")) + ); + DECLARE_LOGGER; +}; +REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictViscoMat_FrictViscoPhys); + +/** 2d functor creating the interaction law (Law2) based on SphereContactGeometry (ScGeom) and FrictViscoPhys of 2 bodies, returning type FrictViscoPM */ +class Law2_ScGeom_FrictViscoPhys_CundallStrackVisco: public LawFunctor{ + public: + OpenMPAccumulator<Real> plasticDissipation; + virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I); + Real elasticEnergy (); + Real getPlasticDissipation(); + void initPlasticDissipation(Real initVal=0); + YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_FrictViscoPhys_CundallStrackVisco,LawFunctor,"Constitutive law for the FrictViscoPM. Corresponds to :yref:`Law2_ScGeom_FrictPhys_CundallStrack` with the only difference that viscous damping in normal direction can be considered.", + ((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)")) + ((bool,sphericalBodies,true,,"If true, compute branch vectors from radii (faster), else use contactPoint-position. Turning this flag true is safe for sphere-sphere contacts and a few other specific cases. It will give wrong values of torques on facets or boxes.")) + ((bool,traceEnergy,false,,"Define the total energy dissipated in plastic slips at all contacts. This will trace only plastic energy in this law, see O.trackEnergy for a more complete energies tracing")) + ((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)")) + ((int,elastPotentialIx,-1,(Attr::hidden|Attr::noSave),"Index for elastic potential energy (with O.trackEnergy)")) + ,, + .def("elasticEnergy",&Law2_ScGeom_FrictViscoPhys_CundallStrackVisco::elasticEnergy,"Compute and return the total elastic energy in all \"FrictViscoPhys\" contacts") + .def("plasticDissipation",&Law2_ScGeom_FrictViscoPhys_CundallStrackVisco::getPlasticDissipation,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:Law2_ScGeom_FrictViscoPhys_CundallStrackVisco::traceEnergy` is true.") + .def("initPlasticDissipation",&Law2_ScGeom_FrictViscoPhys_CundallStrackVisco::initPlasticDissipation,"Initialize cummulated plastic dissipation to a value (0 by default).") + ); + FUNCTOR2D(ScGeom,FrictViscoPhys); + DECLARE_LOGGER; +}; +REGISTER_SERIALIZABLE(Law2_ScGeom_FrictViscoPhys_CundallStrackVisco); + +
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
