# Dissipation power during contact with linear spring dashpot system

Dear all,

I am trying to compute the dissipated power of a granular system in a YADE simulation. In my simulations, I use a linear spring dashpot system (Law2_ScGeom_

Fn = kn delta_n + cn d(delta_n)/dt (normal direction of contact)

Ft = -min(ks delta_t, muFn) (tangential direction of contact)

In the normal direction, only the damper dissipates energy so the normal dissipated power at contact is

Pn = - cn d(delta_n)/dt * d(delta_n)/dt

Considering a collision of two particles with only a normal component (I ran a yade simulation), I have been able to verify that the integrate of Pn during the contact time indeed corresponds to the theoretical energy dissipated for the corresponding restitution coefficient.

My question concerns more the tangential dissipation. Following the same reasoning, I would say that:

if Ft = -muFn:

Pt = -muFn*shearVel

where shearVel is the tangential component of the relative velocity between the two particles. However, in some of the contact laws (those for which the energy dissipation is already coded, for example in ElasticContactL

if Ft = -muFn:

Pt = - muFn * ((Ft[t-1]-ks delta_t) - Ft)/(ks*dt)

where Ft[t-1] is the shear force at previous time step and delta_t corresponds to a the shear increment (I assume delta_t = shearVel*dt). You can find the corresponding portion of code of ElasticContactL

I don't really understand why the sliding velocity has to be computed on the difference of forces and what this physically means. If someone understands that or have a reference to share which gives an explanation I would be grateful.

Thanks in advance,

Rémi

Vector3r& shearForce = geom->rotate(

const Vector3r& shearDisp = geom->shearIncr

shearForce -= phys->ks * shearDisp;

Real maxFs = phys->normalFor

if (!scene-

// PFC3d SlipModel, is using friction angle. CoulombCriterion

if (shearForce.

Real ratio = sqrt(maxFs) / shearForce.norm();

shearForce *= ratio;

}

} else {

//almost the same with additional Vector3r instatinated for energy tracing,

//duplicated block to make sure there is no cost for the instanciation of the vector when traceEnergy==false

if (shearForce.

Real ratio = sqrt(maxFs) / shearForce.norm();

Vector3r trialForce = shearForce; //store prev force for definition of plastic slip

//define the plastic work input and increment the total plastic energy dissipated

shearForce *= ratio;

Real dissip = ((1 / phys->ks) * (trialForce - shearForce)) /*plastic disp*/.

}

}

## Question information

- Language:
- English Edit question

- Status:
- Answered

- For:
- Yade Edit question

- Assignee:
- No assignee Edit question

- Last query:
- 2021-02-10

- Last reply:
- 2021-02-12

Jérôme Duriez (jduriez) said : | #1 |

Hi,

I did not look so closely to your code in your text but, for what concerns the line at [*] the "(1 / phys->ks) * (trialForce - shearForce)" term gives the plastic/

(Not the total one, in an elasto-plastic wording)

It is then just fair to multiply it with current force to get a purely dissipated work.

(and, in case that was your question, there is no need to speak about dashpot ;-) )

## Can you help with this problem?

Provide an answer of your own, or ask Rémi Chassagne for more information if necessary.