how to calculate the plastic part of shear displacement in law2-cohesionMoment

Asked by Huang peilun

Hi,
I want to calculate the plastic shear displacement of any specific contact during every time step. I'm using Law2_ScGeom6D_CohFrictPhys_CohesionMoment[1]. In the source code, the plastic displacement is calculated by the following code (L140-L159):

Vector3r& shearForce = geom->rotate(phys->shearForce);
  const Vector3r& dus = geom->shearIncrement();

  //Linear elasticity giving "trial" shear force
  shearForce -= phys->ks * dus;

  Real Fs = phys->shearForce.norm();
  Real maxFs = phys->shearAdhesion;
  if (!phys->cohesionDisablesFriction || maxFs == 0) maxFs += Fn * phys->tangensOfFrictionAngle;
  maxFs = math::max((Real)0, maxFs);
  if (Fs > maxFs) { //Plasticity condition on shear force
   if (phys->fragile && !phys->cohesionBroken) {
    phys->SetBreakingState();
    maxFs = max((Real)0, Fn * phys->tangensOfFrictionAngle);
   }
   maxFs = maxFs / Fs;
   Vector3r trialForce = shearForce;
   shearForce *= maxFs;
   if (scene->trackEnergy || traceEnergy) {
    Real sheardissip = ((1 / phys->ks) * (trialForce - shearForce)) /*plastic disp*/.dot(shearForce) /*active force*/;

Briefly, the plastic shear displacement is calculated by 1/ks*(trialForce - shearForce). However, it seems trialForce is not available during the simulation process. Though shearInc is available, which is dus in the source code, it is not exactly the same as 1/ks*(trialForce - shearForce). Therefore I want to ask if there is a simple way to get the plastic shear displacement instead of modify the source code?

Thank you!

[1] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/dem/CohesiveFrictionalContactLaw.cpp

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Bruno Chareyre
Solved:
Last query:
Last reply:
Revision history for this message
Best Bruno Chareyre (bruno-chareyre) said :
#1

Hi,
No, I don't see a simple way to track this "per contact".
Energy tracking would give a total for all contacts, not what you want.
You'll have to modify the source code a little.
For instance, add a "Real dissipation" attribute to interaction physics and accumulate there.
I hope it helps.
Bruno

Revision history for this message
Huang peilun (hpl16) said :
#2

Thanks Bruno Chareyre, that solved my question.