A New Constitutive Law: Viscoelastic Burger's model

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();

Hi Behzad,

you can also ask compiler to correct you :-) it would warn you about syntax
errors (ending function } in Implementation), name errors (it would
complain about CohesiveBurgersContactLaw) etc.

Law2::normElastEnergy, Law2::normalDisplacement() etc. are not defined..

Nothing is done in void Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys::go,
you just declare some Reals, and do not use it anywhere (the compiler would
warn you about unused variables)

void Ip2_FrictMat_CohBurgersMat_CohFrictPhys::go calls itself (again
compiler would complain I would guess). See [1] how it is done for
FrictMat_CpmMat.

Law2 is called only once and only values valid globally for the law should
be declared there. Local values (different for each interaction) should go
to IPhys. I feel it is the case of all BurgXX variables, if they depends on
particle size for example. What is the meaning of them in your Law2? aren't
they the same as IPhys:BurXX, which should be computed from material?

IPhys should also store internal variables of the physics if you need some
(I fell like you should use some for the dampers, but maybe it is not
necessary).

So, give the compiler chance to teach you first :-) of course, if you have
a problem, you can ask, but to check everything manually, it is a bit
complex..

cheers
Jan

[1]
http://bazaar.launchpad.net/~yade-pkg/yade/git-trunk/view/head:/pkg/dem/ConcretePM.cpp#L26

2014-10-15 22:46 GMT+02:00 behzad <email address hidden>:

> New question #255780 on Yade:
> https://answers.launchpad.net/yade/+question/255780
>
> Hi there,
>
> It's a while I'm working to derive a new constitutive law which gives the
> properties of viscoelastic Burger's model.
> I share here the first drafts of the model which is integrated in one
> header and one implementation file.
> It's a draft and I have not yet tested it. It's shared here to find the
> possible errors.
> Any comments and helps are highly appreciated.
>
> Thanks
>
>
> HEADER:
>
>
>
>
>
> #pragma once
>
> #include<yade/core/GlobalEngine.hpp>
> #include<yade/pkg/common/Dispatching.hpp>
> #include<yade/pkg/dem/ScGeom.hpp>
> #include<boost/tuple/tuple.hpp>
> #include<yade/pkg/dem/CohesiveFrictionalContactLaw.hpp>
> #include<yade/pkg/dem/FrictPhys.hpp>
>
>
>
> //********************** CohBurgersMat ****************************/
>
> class CohBurgersMat : public CohFrictMat
> {
> public :
> virtual ~CohBurgersMat ();
>
> /// Serialization
> YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohBurgersMat,FrictMat,"",
> ((Real,Em,1e3,,"Stiffness of Maxwell's spring"))
> ((Real,Ek,1e3,,"Stiffness of Kelvin's spring"))
> ((Real,Cm,10,,"Viscosity of Maxwell's dashpot"))
> ((Real,Ck,10,,"Viscosity of Kelvin's dashpot"))
> ((Real,poissonRatio,0.5,,""))
> createIndex();
> );
> /// Indexable
> REGISTER_CLASS_INDEX(CohBurgersMat,CohFrictMat);
> };
>
> REGISTER_SERIALIZABLE(CohBurgersMat);
>
>
>
> //********************** CohBurgersPhys ****************************/
>
> class CohBurgersPhys : public CohFrictPhys
> {
> public :
> virtual ~CohBurgerPhys() {};
> void SetBreakingState() {cohesionBroken = true;
> normalAdhesion = 0; shearAdhesion = 0;};
>
> YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohBurgerPhys,CohFrictPhys,"",
> ((bool,cohesionDisablesFriction,false,,"is shear strength
> the sum of friction and adhesion or only adhesion?"))
> ((bool,cohesionBroken,true,,"is cohesion active? Set to
> false at the creation of a cohesive contact, and set to true when a fragile
> contact is broken"))
> ((bool,fragile,true,,"do cohesion disapear when contact
> strength is exceeded?"))
> ((Real,BurEmn,NaN,,"Stiffness of Maxwell's
> spring(normal)"))
> ((Real,BurEkn,NaN,,"Stiffness of Kelvin's spring(normal)"))
> ((Real,BurCmn,NaN,,"Viscosity of Maxwell's
> dashpot(normal)"))
> ((Real,BurCkn,NaN,,"Viscosity of Kelvin's
> dashpot(normal)"))
> ((Real,BurEms,NaN,,"Stiffness of Maxwell's spring(shear)"))
> ((Real,BurEks,NaN,,"Stiffness of Kelvin's spring(shear)"))
> ((Real,BurCms,NaN,,"Viscosity of Maxwell's
> dashpot(shear)"))
> ((Real,BurCks,NaN,,"Viscosity of Kelvin's dashpot(shear)"))
> ((Real,poissonRatio,NaN,,""))
> ((Real,normalAdhesion,0,,"tensile strength"))
> ((Real,shearAdhesion,0,,"cohesive part of the shear
> strength (a frictional term might be added depending on
> :yref:`CohFrictPhys::cohesionDisablesFriction`)"))
> ,
> createIndex();
> );
> /// Indexable
> REGISTER_CLASS_INDEX(CohBurgerPhys,CohFrictPhys);
> };
>
> REGISTER_SERIALIZABLE(CohBurgerPhys);
>
>
>
>
> //********************** 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,
>
> ((bool,setCohesionNow,false,,"If true, assign cohesion to
> all existing contacts in current time-step. The flag is turned false
> automatically, so that assignment is done in the current timestep only."))
> ((bool,setCohesionOnNewContacts,false,,"If true, assign
> cohesion at all new contacts. If false, only existing contacts can be
> cohesive (also see
> :yref:`Ip2_CohFrictMat_CohFrictMat_CohFrictPhys::setCohesionNow`), and new
> contacts are only frictional."))
> ,
> cohesionDefinitionIteration = -1;
> );
> FUNCTOR2D(CohBurgersMat,CohBurgersMat);
> };
>
> REGISTER_SERIALIZABLE(Ip2_CohBurgersMat_CohBurgersMat_CohBurgersPhys);
>
>
>
> //********************** Ip2_FrictMat_CohBurgersMat_CohFrictPhys
> ****************************/
>
>
> class Ip2_FrictMat_CohBurgersMat_CohFrictPhys : 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_CohFrictPhys,IPhysFunctor,
>
> ((bool,setCohesionNow,false,,"If true, assign cohesion to
> all existing contacts in current time-step. The flag is turned false
> automatically, so that assignment is done in the current timestep only."))
> ((bool,setCohesionOnNewContacts,false,,"If true, assign
> cohesion at all new contacts. If false, only existing contacts can be
> cohesive (also see
> :yref:`Ip2_CohFrictMat_CohFrictMat_CohFrictPhys::setCohesionNow`), and new
> contacts are only frictional."))
> ,
> cohesionDefinitionIteration = -1;
> );
> FUNCTOR2D(FrictMat,CohBurgersMat);
> };
>
> REGISTER_SERIALIZABLE(Ip2_FrictMat_CohBurgersMat_CohFrictPhys);
>
>
>
>
> //******************** Law2_ScGeom_CohBurgersPhys_CohesiveBurgers
> *************************/
>
>
> class Law2_ScGeom_CohBurgersPhys_CohesiveBurgers: public LawFunctor{
> public:
> OpenMPAccumulator<Real> plasticDissipation;
> Real normElastEnergy();
> Real shearElastEnergy();
> Real normalDisplacement();
> Real shearDisplacement();
> Real totalElastEnergy();
> Real ((BurgEmn, 0.0, "Em in normal direction"))
> Real ((BurgEkn, 0.0, "Ek in normal direction"))
> Real ((BurgEms, 0.0, "Em in shear direction"))
> Real ((BurgEks, 0.0, "Ek in shear direction"))
> Real ((BurgCmn, 0.0, "Cm in normal direction"))
> Real ((BurgCkn, 0.0, "Ck in normal direction"))
> Real ((BurgCms, 0.0, "Cm in shear direction"))
> Real ((BurgCks, 0.0, "Ck in shear direction"))
> Real ((BurgPoi,0.0,"Poisson ratio"))
>
> void initPlasticDissipation(Real initVal=0);
> virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>&
> _phys, Interaction* I);
>
> YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_CohBurgerPhys_CohesiveBurger,
> ((bool,neverErase,false,,))
> ,,
>
> .def("normElastEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::normElastEnergy,"Compute
> normal elastic energy.")
>
> .def("shearElastEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::shearElastEnergy,"Compute
> shear elastic energy.")
>
> .def("normalDisplacement",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::normalDisplacement,"Compute
> normal displacement.")
>
> .def("shearDisplacement",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::shearDisplacement,"Compute
> shear displacement.")
>
> .def("elasticEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::totalElastEnergy,"Compute
> total elastic energy.")
> );
> FUNCTOR2D(ScGeom,CohBurgersPhys);
> DECLARE_LOGGER;
> };
> REGISTER_SERIALIZABLE(Law2_ScGeom_CohBurgersPhys_CohesiveBurgers);
>
>
> --
> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hi,

I saved the new files(burgersmodl.hpp and burgersmodel.cpp) in trunk/pkg/dem.
Compilation is not going through. I guess I need to modify something in make files, right?

Hi, did you run cmake before make?
Jan

2014-10-17 17:16 GMT+02:00 behzad <email address hidden>:

> Question #255780 on Yade changed:
> https://answers.launchpad.net/yade/+question/255780
>
> Status: Answered => Open
>
> behzad is still having a problem:
>
> Hi,
>
> I saved the new files(burgersmodl.hpp and burgersmodel.cpp) in
> trunk/pkg/dem.
> Compilation is not going through. I guess I need to modify something in
> make files, right?
>
> --
> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

no!

Can I just compile the new files and not all the source code again?

Since your files are surely deeply related with (i.e. include) lot of other source files, I guess it is not possible.

But new compilation of yade in a whole should involve only what has changed : so your new files and all other files that include them. (though, not 100% sure I'm right in your case, when you run a new cmake)

Jerome

Hi Behzad,
just now I did some changes in source code, run cmake and than make and
only the changes are compiled (of course not just the new files, but also
files and actions dependent on them). IMO the same case should be for you..
cheers
Jan

2014-10-17 17:41 GMT+02:00 behzad <email address hidden>:

> Question #255780 on Yade changed:
> https://answers.launchpad.net/yade/+question/255780
>
> behzad posted a new comment:
> no!
>
> Can I just compile the new files and not all the source code again?
>
> --
> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

I ran

cmake
make

then I restarted yade. But when I asked for the list of materials as:

yade.system.childClasses('Material')

the new class of material which I defined in the new files is not there. So, the changes are not made in the code.

make install

2014-10-17 18:21 GMT+02:00 behzad <email address hidden>:

> Question #255780 on Yade changed:
> https://answers.launchpad.net/yade/+question/255780
>
> behzad posted a new comment:
> I ran
>
> cmake
> make
>
> then I restarted yade. But when I asked for the list of materials as:
>
> yade.system.childClasses('Material')
>
>
> the new class of material which I defined in the new files is not there.
> So, the changes are not made in the code.
>
> --
> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Thank you!