Hi there, I have developed a new contact model, namely Burger's viscoelastic model for YADE. I have tried and verified the model in my own YADE trunk. I'd like to publish it on the main trunk to get it reviewed and further developed (if any)
was not able to commit the changes directly on github, using git push command. My username is bmajidi However, I have attached the model files to this email and I appreciate if you can commit the changes for me. Thank you, -- Behzad Majidi Research Student Aluminium Research Centre REGAL-Université Laval Québec, QC, Canada Please consider the environment before printing this E-mail. -- Behzad Majidi Research Student Aluminium Research Centre REGAL-Université Laval Québec, QC, Canada Please consider the environment before printing this E-mail.
// 2014 © Behzad Majidi <[email protected]> // 2014 © Ricardo Pieralisi <[email protected]> #include"Cohburgers.hpp" #include<core/State.hpp> #include<core/Interaction.hpp> #include<pkg/dem/ScGeom.hpp> #include<core/Omega.hpp> #include<core/Scene.hpp> #include<pkg/dem/Shop.hpp> #include<pkg/common/InteractionLoop.hpp> #include<pkg/common/Sphere.hpp> #include<pkg/common/Facet.hpp> #include<pkg/common/Wall.hpp> YADE_PLUGIN((CohBurgersMat)(CohBurgersPhys)(Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys)(Ip2_CohFrictMat_CohBurgersMat_CohBurgersPhys)(Ip2_FrictMat_CohBurgersMat_CohBurgersPhys)(Law2_ScGeom6D_CohBurgersPhys_CohesiveBurgers)); //********************** Ip2_FrictMat_CohBurgersMat_CohBurgersPhys ****************************/ void Ip2_FrictMat_CohBurgersMat_CohBurgersPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) { if(interaction->phys) return; shared_ptr<CohBurgersPhys> phys (new CohBurgersPhys()); int i1 = b1->getClassIndex(); // get actual material index of b1 int i2 = b2->getClassIndex(); // get actual material index of b2 int cbmi = CohBurgersMat::getClassIndexStatic(); // get index of CohBurgersMat FrictMat* mat1; CohBurgersMat* mat2; if (i2 == cbmi) { // b2 is CohburgersMat and b1 is FrictMat mat1 = dynamic_cast<FrictMat*>(b1.get()); mat2 = dynamic_cast<CohBurgersMat*>(b2.get()); } else if (i1 == cbmi) { // b1 is CohburgersMat and b2 is FrictMat mat1 = dynamic_cast<FrictMat*>(b2.get()); mat2 = dynamic_cast<CohBurgersMat*>(b1.get()); } else { // should not happen, but to be sure.. LOG_FATAL("TODO"); } ScGeom6D* geom= YADE_CAST<ScGeom6D*>(interaction->geom.get()); Real D1 = geom->radius1*2; Real D2 = geom->radius2*2; //Load the rheological parameters Real kn1 = (mat1->young)*D1; Real ks1 = (mat1->young)*D1*(mat1->poisson); Real kmn2 = (mat2->kmn)*D2; Real kkn2 = (mat2->kkn)*D2; Real cmn2 = (mat2->cmn)*D2; Real ckn2 = (mat2->ckn)*D2; Real kms2 = (mat2->kms)*D2; Real kks2 = (mat2->kks)*D2; Real cms2 = (mat2->cms)*D2; Real cks2 = (mat2->cks)*D2; phys->kmn = kn1; phys->kms = ks1; phys->cmn = cmn2; phys->cms = cms2; phys->kkn = kkn2; phys->kks = kks2; phys->ckn = ckn2; phys->cks = cks2; if ((setCohesionOnNewContacts || setCohesionNow) && mat2->isCohesive) { phys->MaxDef = 0.125*(D1+D2); phys->normalAdhesion= (mat2->normalCohesion)*0.5*pow(std::min(D1*0.5,D2*0.5),2); phys->shearAdhesion= (mat2->shearCohesion)*0.5*pow(std::min(D1*0.5,D2*0.5),2); phys->fragile=mat2->fragile; phys->cohesionBroken = false; phys->initCohesion = true; } interaction->phys = shared_ptr<CohBurgersPhys>(phys); } //********************** Ip2_CohFrictMat_CohBurgersMat_CohBurgersPhys ****************************/ void Ip2_CohFrictMat_CohBurgersMat_CohBurgersPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) { if(interaction->phys) return; shared_ptr<CohBurgersPhys> phys (new CohBurgersPhys()); int i1 = b1->getClassIndex(); // get actual material index of b1 int i2 = b2->getClassIndex(); // get actual material index of b2 int cbmi = CohBurgersMat::getClassIndexStatic(); // get index of CohBurgersMat CohFrictMat* mat1; CohBurgersMat* mat2; if (i2 == cbmi) { // b2 is CohburgersMat and b1 is FrictMat mat1 = dynamic_cast<CohFrictMat*>(b1.get()); mat2 = dynamic_cast<CohBurgersMat*>(b2.get()); } else if (i1 == cbmi) { // b1 is CohburgersMat and b2 is FrictMat mat1 = dynamic_cast<CohFrictMat*>(b2.get()); mat2 = dynamic_cast<CohBurgersMat*>(b1.get()); } else { // should not happen, but to be sure.. LOG_FATAL("TODO"); } ScGeom6D* geom= YADE_CAST<ScGeom6D*>(interaction->geom.get()); Real D1 = geom->radius1*2; Real D2 = geom->radius2*2; //Load the rheological parameters Real kn1 = (mat1->young)*D1; Real ks1 = (mat1->young)*D1*(mat1->poisson); Real kmn2 = (mat2->kmn)*D2; Real kkn2 = (mat2->kkn)*D2; Real cmn2 = (mat2->cmn)*D2; Real ckn2 = (mat2->ckn)*D2; Real kms2 = (mat2->kms)*D2; Real kks2 = (mat2->kks)*D2; Real cms2 = (mat2->cms)*D2; Real cks2 = (mat2->cks)*D2; phys->kmn = kn1; phys->kms = ks1; phys->cmn = cmn2; phys->cms = cms2; phys->kkn = kkn2; phys->kks = kks2; phys->ckn = ckn2; phys->cks = cks2; if ((setCohesionOnNewContacts || setCohesionNow) && mat1->isCohesive && mat2->isCohesive) { phys->MaxDef = 0.125*(D1+D2); phys->normalAdhesion= (mat2->normalCohesion)*0.5*pow(std::min(D1*0.5,D2*0.5),2); phys->shearAdhesion= (mat2->shearCohesion)*0.5*pow(std::min(D1*0.5,D2*0.5),2); phys->fragile=mat2->fragile; phys->cohesionBroken = false; phys->initCohesion = true; } interaction->phys = shared_ptr<CohBurgersPhys>(phys); } //********************** Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys ****************************/ void Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) { if(interaction->phys) return; shared_ptr<CohBurgersPhys> phys (new CohBurgersPhys()); CohBurgersMat* mat1 = YADE_CAST<CohBurgersMat*>(b1.get()); CohBurgersMat* mat2 = YADE_CAST<CohBurgersMat*>(b2.get()); ScGeom6D* geom= YADE_CAST<ScGeom6D*>(interaction->geom.get()); Real D1 = geom->radius1*2; Real D2 = geom->radius2*2; //Load the rheological parameters Real kmn1 = (mat1->kmn)*2*D1; Real kkn1 = (mat1->kkn)*2*D1; Real cmn1 = (mat1->cmn)*2*D1; Real ckn1 = (mat1->ckn)*2*D1; Real kms1 = (mat1->kms)*2*D1; Real kks1 = (mat1->kks)*2*D1; Real cms1 = (mat1->cms)*2*D1; Real cks1 = (mat1->cks)*2*D1; Real kmn2 = (mat2->kmn)*2*D2; Real kkn2 = (mat2->kkn)*2*D2; Real cmn2 = (mat2->cmn)*2*D2; Real ckn2 = (mat2->ckn)*2*D2; Real kms2 = (mat2->kms)*2*D2; Real kks2 = (mat2->kks)*2*D2; Real cms2 = (mat2->cms)*2*D2; Real cks2 = (mat2->cks)*2*D2; phys->kmn = contactParameterCalculation(kmn1,kmn2); phys->kms = contactParameterCalculation(kms1,kms2); phys->cmn = contactParameterCalculation(cmn1,cmn2); phys->cms = contactParameterCalculation(cms1,cms2); phys->kkn = contactParameterCalculation(kkn1,kkn2); phys->kks = contactParameterCalculation(kks1,kks2); phys->ckn = contactParameterCalculation(ckn1,ckn2); phys->cks = contactParameterCalculation(cks1,cks2); if ((setCohesionOnNewContacts || setCohesionNow) && mat1->isCohesive && mat2->isCohesive) { phys->MaxDef = 0.25*(D1+D2); phys->normalAdhesion= std::min(mat1->normalCohesion, mat2->normalCohesion)*pow(std::min(D1*0.5,D2*0.5),2); phys->shearAdhesion= (mat2->shearCohesion)*0.5*pow(std::min(D1*0.5,D2*0.5),2); phys->fragile=(mat1->fragile || mat2->fragile); phys->cohesionBroken = false; phys->initCohesion = true; } interaction->phys = shared_ptr<CohBurgersPhys>(phys); } //********************** CohesiveBurgersContactLaw ****************************/ bool Law2_ScGeom6D_CohBurgersPhys_CohesiveBurgers::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I) { Vector3r force = Vector3r::Zero(); Vector3r torque1 = Vector3r::Zero(); Vector3r torque2 = Vector3r::Zero(); CohBurgersPhys& phys=*YADE_CAST<CohBurgersPhys*>(_phys.get()); if (computeForceTorqueCohBurgers(_geom, _phys, I, force, torque1, torque2) and (I->isActive)){ const int id1 = I->getId1(); const int id2 = I->getId2(); addForce (id1,-force,scene); addForce (id2, force,scene); addTorque(id1, torque1,scene); addTorque(id2, torque2,scene); return true; } else return false; } bool computeForceTorqueCohBurgers(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r & torque2) { const ScGeom6D& geom=*YADE_CAST<ScGeom6D*>(_geom.get()); CohBurgersPhys& phys=*YADE_CAST<CohBurgersPhys*>(_phys.get()); Scene* scene=Omega::instance().getScene().get(); const int id1 = I->getId1(); const int id2 = I->getId2(); const BodyContainer& bodies = *scene->bodies; const State& de1 = *YADE_CAST<State*>(bodies[id1]->state.get()); const State& de2 = *YADE_CAST<State*>(bodies[id2]->state.get()); const Real& dt = scene->dt; if (phys.fragile && phys.normalAdhesion > 0 && phys.normalForce.norm() > phys.normalAdhesion) { phys.cohesionBroken= true; return false; } if (phys.fragile && std::abs(geom.penetrationDepth) > phys.MaxDef) { phys.cohesionBroken= true; return false; } if (I->isFresh(scene)) { phys.ukn = Vector3r(0,0,0); phys.uks = Vector3r(0,0,0); phys.shearForce = Vector3r(0,0,0); phys.normalForce = phys.kmn * geom.penetrationDepth * geom.normal; return true; } else { Real An = 1.0 + phys.kkn * dt / (2.0 * phys.ckn); Real Bn = 1.0 - phys.kkn * dt / (2.0 * phys.ckn); Real Cn = dt / (2.0 * phys.ckn * An) + 1.0 / phys.kmn + dt / (2.0 * phys.cmn); Real Dn = dt / (2.0 * phys.ckn * An) - 1.0 / phys.kmn + dt / (2.0 * phys.cmn); Real As = 1.0 + phys.kks * dt / (2.0 * phys.cks); Real Bs = 1.0 - phys.kks * dt / (2.0 * phys.cks); Real Cs = dt / (2.0 * phys.cks * As) + 1.0 / phys.kms + dt / (2.0 * phys.cms); Real Ds = dt / (2.0 * phys.cks * As) - 1.0 / phys.kms + dt / (2.0 * phys.cms); const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero(); const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(I->cellDist): Vector3r::Zero(); const Vector3r c1x = (geom.contactPoint - de1.pos); const Vector3r c2x = (geom.contactPoint - de2.pos - shift2); const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) + shiftVel; const Vector3r normalVelocity = geom.normal.dot(relativeVelocity) * geom.normal; const Vector3r shearVelocity = relativeVelocity-normalVelocity; const Vector3r normalForceT = (normalVelocity * dt + phys.ukn * (1.0 - Bn / An) - phys.normalForce * Dn) / Cn; phys.ukn = (phys.ukn * Bn + (normalForceT + phys.normalForce) * dt / (2.0 * phys.ckn)) / An; phys.normalForce = normalForceT; const Vector3r shearForceT = -1 * (shearVelocity * dt + phys.uks * (1.0 - Bs / As) - phys.shearForce * Ds) / Cs; phys.uks = (phys.uks * Bs + (shearForceT + phys.shearForce) * dt / (2.0 * phys.cks)) / As; phys.shearForce = shearForceT; const Real maxFs = phys.normalForce.squaredNorm() * std::pow(phys.tangensOfFrictionAngle,2); if( phys.shearForce.squaredNorm() > maxFs ) { const Real ratio = sqrt(maxFs) / phys.shearForce.norm(); phys.shearForce *= ratio; } force = phys.normalForce + phys.shearForce; torque1 = -c1x.cross(force); torque2 = c2x.cross(force); return true; } }
// 2014 © Behzad Majidi <[email protected]> // 2014 © Ricardo Pieralisi <[email protected]> #pragma once #include<core/Material.hpp> #include<pkg/dem/FrictPhys.hpp> #include<pkg/dem/ViscoelasticPM.hpp> #include<pkg/dem/CohesiveFrictionalContactLaw.hpp> #include<pkg/common/Dispatching.hpp> #include<pkg/dem/ScGeom.hpp> #include<pkg/dem/DemXDofGeom.hpp> #include<pkg/common/MatchMaker.hpp> //********************** CohBurgersMat ****************************/ class CohBurgersMat : public CohFrictMat { public : virtual ~CohBurgersMat () {}; /// Serialization YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohBurgersMat,CohFrictMat,"", ((Real,kmn,NaN,,"Normal elastic stiffness for Maxwell.")) ((Real,cmn,NaN,,"Normal viscous constant for Maxwell.")) ((Real,kms,NaN,,"Shear elastic stiffness for Maxwell.")) ((Real,cms,NaN,,"Shear viscous constant for Maxwell.")) ((Real,kkn,NaN,,"Normal elastic stiffness for Kelvin.")) ((Real,ckn,NaN,,"Normal viscous constant for Kelvin.")) ((Real,kks,NaN,,"Shear elastic stiffness for Kelvin.")) ((Real,cks,NaN,,"Shear viscous constant for Kelvin.")), createIndex(); ); /// Indexable REGISTER_CLASS_INDEX(CohBurgersMat,CohFrictMat); }; REGISTER_SERIALIZABLE(CohBurgersMat); //********************** CohBurgersPhys ****************************/ class CohBurgersPhys : public CohFrictPhys { public : virtual ~CohBurgersPhys() {}; void SetBreakingState() {cohesionBroken = true; normalAdhesion = 0; shearAdhesion = 0;}; YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohBurgersPhys,CohFrictPhys,"", ((Real,kmn,NaN,,"Stiffness of Maxwell's spring(normal)")) ((Real,kkn,NaN,,"Stiffness of Kelvin's spring(normal)")) ((Real,cmn,NaN,,"Viscosity of Maxwell's dashpot(normal)")) ((Real,ckn,NaN,,"Viscosity of Kelvin's dashpot(normal)")) ((Real,kms,NaN,,"Stiffness of Maxwell's spring(shear)")) ((Real,kks,NaN,,"Stiffness of Kelvin's spring(shear)")) ((Real,cms,NaN,,"Viscosity of Maxwell's dashpot(shear)")) ((Real,cks,NaN,,"Viscosity of Kelvin's dashpot(shear)")) ((Real,MaxDef,NaN,,"maximum permitted deformation")) ((Vector3r,ukn,Vector3r(0,0,0),,"Normal displacement")) ((Vector3r,uks,Vector3r(0,0,0),,"Shear displacement")), createIndex(); ) REGISTER_CLASS_INDEX(CohBurgersPhys,CohFrictPhys); }; REGISTER_SERIALIZABLE(CohBurgersPhys); //********************** Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys ****************************/ class Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys: public IPhysFunctor { public : virtual void go( const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction); int cohesionDefinitionIteration; YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys,IPhysFunctor, "Convert 2 instances of :yref:`CohBurgersMat` to :yref:`CohBurgersPhys` using the rule of consecutive connection.", ((bool,setCohesionNow,false,,"If true, assign cohesion to all existing contacts")) ((bool,setCohesionOnNewContacts,false,,"If true, assign cohesion at all new contacts.")) , cohesionDefinitionIteration = -1; ); FUNCTOR2D(CohBurgersMat,CohBurgersMat); }; REGISTER_SERIALIZABLE(Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys); class Ip2_CohFrictMat_CohBurgersMat_CohBurgersPhys: public IPhysFunctor { public : virtual void go( const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction); int cohesionDefinitionIteration; YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_CohFrictMat_CohBurgersMat_CohBurgersPhys,IPhysFunctor, "converts an interaction of FrictMat and CohBurgersMat to CohBurgersPhys", ((bool,setCohesionNow,false,,"If true, assign cohesion to all existing contacts")) ((bool,setCohesionOnNewContacts,false,,"If true, assign cohesion at all new contacts.")) , cohesionDefinitionIteration = -1; ); FUNCTOR2D(CohFrictMat,CohBurgersMat); }; REGISTER_SERIALIZABLE(Ip2_CohFrictMat_CohBurgersMat_CohBurgersPhys); class Ip2_FrictMat_CohBurgersMat_CohBurgersPhys: public IPhysFunctor { public : virtual void go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction); int cohesionDefinitionIteration; YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_FrictMat_CohBurgersMat_CohBurgersPhys,IPhysFunctor, "converts an interaction of FrictMat and CohBurgersMat to CohBurgersPhys", ((bool,setCohesionNow,false,,"If true, assign cohesion to all existing contacts")) ((bool,setCohesionOnNewContacts,false,,"If true, assign cohesion at all new contacts.")) , cohesionDefinitionIteration = -1; ); FUNCTOR2D(FrictMat,CohBurgersMat); }; REGISTER_SERIALIZABLE(Ip2_FrictMat_CohBurgersMat_CohBurgersPhys); //******************** Law2_ScGeom6D_CohBurgersPhys_CohesiveBurgers *************************/ class Law2_ScGeom6D_CohBurgersPhys_CohesiveBurgers: public LawFunctor { public: virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction*); FUNCTOR2D(ScGeom6D,CohBurgersPhys); YADE_CLASS_BASE_DOC(Law2_ScGeom6D_CohBurgersPhys_CohesiveBurgers, LawFunctor,""); }; REGISTER_SERIALIZABLE(Law2_ScGeom6D_CohBurgersPhys_CohesiveBurgers); bool computeForceTorqueCohBurgers(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r & torque2);
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
