Info on traceEnergy and O.trackEnergy

Asked by Vasileios Angelidakis on 2019-10-06

Hi,

I am currently developing a "traceEnergy" feature in the Law2_SCG_KnKsPBPhys_KnKsPBLaw for the Potential Blocks and I could use some advice. :) I am trying to make the implementation as close to existing ones for other constitutive laws, as possible, in order to keep programming similarity and make future maintenance easier, so I would like to have your opinion.

#Question 1:
To calculate the plastic dissipation due to friction for a sliding contact, I see in the Law2_ScGeom_FrictViscoPhys_CundallStrackVisco [1]:

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

The same exists in other constitutive laws, e.g. in ElasticContactLaw.cpp and Polyhedra.cpp.
Why do we only check the sign of "dissip" when we are about to add it in the O.energy module (i.e. traceEnergy=False), but not when we are about to add it in the OpenMPAccumulator "plasticDissipation" (i.e. traceEnergy=true)?

#Question 2:
In the Law2_SCG_KnKsPBPhys_KnKsPBLaw, we consider dissipation due to viscous damping in the normal direction. Would you suggest I create an OpemMPAccumulator "normDampDissip", like in the HertzMindlin implementation [2] and use "traceEnergy" to decide where to store this, like above (i.e. either in the OpemMPAccumulator or in O.energy), or was "traceEnergy" conceived to be used only to decide where to store the plastic dissipation due to friction?

Also, when I add this dissipation to O.energy, should I check its sign, like we do for "dissip" above, and do you think it would be more convenient to reset the value in every timestep, or monitor the cumulative value (i.e. reset=false), like we do for the dissipation due to friction?

#Question 3:
In the Law2_SCG_KnKsPBPhys_KnKsPBLaw, after we calculate the normalViscous force, we deduct it from the normal force [3]:

 phys->normalForce -= phys->normalViscous;

When I calculate the elastic potential energy, should I consider the reduced potential energy or the originally elastic one? i.e.
Reduced potential energy:
 scene->energy->add(0.5*(phys->normalForce.squaredNorm()/phys->kn+phys->shearForce.squaredNorm()/phys->ks),"elastPotential",elastPotentialIx,/*reset*/true);

Elastic potential energy :
 scene->energy->add(0.5*((phys->normalForce+phys->normalViscous).squaredNorm()/phys->kn+(phys->shearForce+phys->shearViscous).squaredNorm()/phys->ks),"elastPotential",elastPotentialIx,/*reset*/true);

Apologies in advance for not opening three different questions; I just feel these should be in one thread, as they concern a very specific part of the source code.

Many thanks,
Vasileios

[1] https://gitlab.com/yade-dev/trunk/blob/master/pkg/dem/FrictViscoPM.cpp#L210-212
[2] https://gitlab.com/yade-dev/trunk/blob/master/pkg/dem/HertzMindlin.cpp#L365
[3] https://gitlab.com/yade-dev/trunk/blob/master/pkg/dem/KnKsPBLaw.cpp#L228

Question information

Language:
English Edit question
Status:
Open
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2019-10-06
Last reply:

Can you help with this problem?

Provide an answer of your own, or ask Vasileios Angelidakis for more information if necessary.

To post a message you must log in.