Question #255780 on Yade changed:
https://answers.launchpad.net/yade/+question/255780
behzad posted a new comment:
IMPLEMENTATION
#include<yade/pkg/dem/Burgers.hpp>
#include<yade/pkg/dem/ScGeom.hpp>
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
//********************** 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)
{
CohBurgersMat* sdec1 = static_cast<CohBurgersMat*>(b1.get());
CohBurgersMat* sdec2 = static_cast<CohBurgersMat*>(b2.get());
ScGeom* geom = YADE_CAST<ScGeom*>(interaction->geom.get());
if (setCohesionNow && cohesionDefinitionIteration==-1)
cohesionDefinitionIteration=scene->iter;
if (setCohesionNow && cohesionDefinitionIteration!=-1 &&
cohesionDefinitionIteration!=scene->iter) {
cohesionDefinitionIteration = -1;
setCohesionNow = 0;}
if (geom) {
if (!interaction->phys) {
interaction->phys = shared_ptr<CohBurgersPhys>(new
CohBurgersPhys());
CohBurgersPhys* contactPhysics =
YADE_CAST<CohBurgersPhys*>(interaction->phys.get());
Real Da = geom->radius1;
Real Db = geom->radius2;
Real BurEmn =sdec1->Em;
Real BurEkn =sdec1->Ek;
Real BurCmn =sdec1->Cm;
Real BurCkn =sdec1->Ck;
Real BurEms =sdec1->poissonRatio*BurEmn;
Real BurEks =sdec1->poissonRatio*BurEkn;
Real BurCms =sdec1->poissonRatio*BurCmn;
Real BurCks =sdec1->poissonRatio*BurCkn;
Real poissonRatio =sdec1->poissonRatio
Real frictionAngle =sdec1->frictionAngle
Real normalAdhesion =sdec1->normalAdhesion
Real shearAdhesion =sdec1->shearAdhesion
//********************** Ip2_FrictMat_CohBurgersMat_CohFrictPhys
****************************/
CREATE_LOGGER(Ip2_FrictMat_CohBurgersMat_CohFrictPhys);
void Ip2_FrictMat_CohBurgersMat_CohFrictPhys::go(const shared_ptr<Material>&
pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>&
interaction){
TIMING_DELTAS_START();
const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(pp1);
const shared_ptr<CohBurgersMat>& mat2 =
YADE_PTR_CAST<CohBurgersMat>(pp2);
Ip2_FrictMat_CohBurgersMat_CohFrictPhys().go(mat1,mat2,interaction);
TIMING_DELTAS_CHECKPOINT("end of Ip2_CohFritPhys");
//********************** CohesiveBurgersContactLaw ****************************/
YADE_PLUGIN((CohesiveBurgersContactLaw)(Law2_ScGeom_CohBurgersPhys)(CohBurgerstMat)(CohBurgersPhys)(Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys));
CREATE_LOGGER(Law2_ScGeom_CohBurgersPhys);
void CohesiveBurgersContactLaw::action()
bool Law2_ScGeom_CohBurgersPhys::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>&
ip, Interaction* contact){
int id1 = contact->getId1(), id2 = contact->getId2();
ScGeom* geom= static_cast<ScGeom*>(ig.get());
FrictPhys* phys = static_cast<CohBurgersPhys*>(ip.get());
if(geom->penetrationDepth <0){
if (neverErase) {
Real Da = geom->radius1;
Real Db = geom->radius2;
Vector3r& shearForce = Phys->shearForce;
Vector3r& normalForce = Phys->normalForce;
Vector3r& ukkn = Phys->ukkn;
Vector3r& shearForceT1 = Phys->shearForceT1;
Vector3r& normalForceT1 = Phys->normalForceT1;
Vector3r& ukks = Phys->ukks;
else return false;
int& Ite = Phys->Ite;
Real Dinit = Da+Db;
Real knk = BurEkn*Dinit;
Real cnk = BurCkn*Dinit;
Real knm = BurEmn*Dinit;
Real cnm = BurCmn*Dinit;
Real ksk = BurEks*Dinit/(2.0+2.0*BurPoi);
Real csk = BurCks*Dinit/(2.0+2.0*BurPoi);
Real ksm = BurEms*Dinit/(2.0+2.0*BurPoi);
Real csm = BurCms*Dinit/(2.0+2.0*BurPoi);
Real un=geom->penetrationDepth;
Real burnA=1.0+knk*dt/(2.0*cnk);
Real burnB=1.0-knk*dt/(2.0*cnk);
Real burnC=dt/(2.0*cnk*burnA)+1.0/knm+dt/(2.0*cnm);
Real burnD=dt/(2.0*cnk*burnA)-1.0/knm+dt/(2.0*cnm);
Real bursA=1.0+ksk*dt/(2.0*csk);
Real bursB=1.0-ksk*dt/(2.0*csk);
Real bursC=dt/(2.0*csk*bursA)+1.0/ksm+dt/(2.0*csm);
Real bursD=dt/(2.0*csk*bursA)-1.0/ksm+dt/(2.0*csm);
if (contact->isNew)
{ ukkn=Vector3r(0,0,0);
ukks=Vector3r(0,0,0);
shearForce=Vector3r(0,0,0);
normalForce=knm*std::max(un,(Real)
0)*geom->normal;
}
Vector3r axis;
Real angle;
axis = Phys->prevNormal.Cross(geom->normal);
shearForce -= shearForce.Cross(axis);
Vector3r summaryAngularVelocity(0,0,0);
if (isDynamic1) summaryAngularVelocity += de1->angularVelocity;
if (isDynamic2) summaryAngularVelocity += de2->angularVelocity;
angle = dt*0.5*geom->normal.Dot(summaryAngularVelocity);
axis = angle*geom->normal;
shearForce -= shearForce.Cross(axis);
Vector3r x = geom->contactPoint;
Vector3r c1x = (x - de1-se3.position);
Vector3r c2x = (x - de2-se3.position);
Vector3r _c1x_= (isDynamic1) ?
geom->radius1*geom->normal : x -de1->zeroPoint;
Vector3r _c2x_= (isDynamic2) ? -
geom->radius2*geom->normal : x -
de2->zeroPoint;
Vector3r relativeVelocity= (de2-
>velocity+de2->angularVelocity.Cross(_c2x_))-(de1->velocity+de1-
>angularVelocity.Cross(_c1x_));
Vector3r normalVelocity =
std::abs(geom->normal.Dot(relativeVelocity))*geom->normal;
Vector3r shearVelocity =relativeVelocity-normalVelocity;
Vector3r shearDisplacement = shearVelocity*dt;
Vector3r normalDisplacement= normalVelocity*dt;
if (contact->isNew)
Ite=Omega::instance().getCurrentIteration();
normalForceT1 = (normalDisplacement+ukkn*(1.0-
burnB/burnA)-normalForce*burnD)/burnC;
ukkn =
(ukkn*burnB+(normalForceT1+normalForce)*dt/(2.0*cnk))/burnA;
normalForce = normalForceT1;
shearForceT1 = -1.0*(shearDisplacement+ukks*(1.0-
bursB/bursA)+bursD*shearForce)/bursC;
ukks = (ukks*bursB-
(shearForceT1+shearForce)*dt/(2.0*csk))/bursA;
shearForce = shearForceT1;
Real maxFs = normalForce.SquaredLength() *
std::pow(Phys->tangensOfFrictionAngle,2);
if( shearForce.SquaredLength() > maxFs )
{
maxFs = Mathr::Sqrt(maxFs) /
shearForce.Length();
shearForce *= maxFs;
}
Vector3r f = normalForce +
shearForce;
static_cast<Force*> ( ncb->physicalActions-
>find( id1 , actionForceIndex).get() )->force -= f;
static_cast<Force*> ( ncb->physicalActions-
>find( id2 , actionForceIndex ).get() )->force += f;
static_cast<Momentum*>( ncb->physicalActions-
>find( id1 , actionMomentumIndex ).get() )->momentum -= c1x.Cross(f);
static_cast<Momentum*>( ncb->physicalActions-
>find( id2 , actionMomentumIndex ).get() )->momentum += c2x.Cross(f);
Phys->prevNormal =geom->normal;
}
}
}
YADE_PLUGIN();
--
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
