Dissipation power during contact with linear spring dashpot system

Asked by Rémi Chassagne on 2021-02-10

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_ViscElPhys_Basic()), that is to say:
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 ElasticContactLaw.cpp), the tangential dissipation is computed differently. It is not the shearVel velocity which is considered but a velocity based on the excess of elastic force compared to the sliding force:
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 ElasticContactLaw.cpp just below.

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.

Rémi

Vector3r& shearForce = geom->rotate(phys->shearForce);
const Vector3r& shearDisp = geom->shearIncrement();
shearForce -= phys->ks * shearDisp;
Real maxFs = phys->normalForce.squaredNorm() * math::pow(phys->tangensOfFrictionAngle, 2);

if (!scene->trackEnergy && !traceEnergy) { //Update force but don't compute energy terms (see below))
// PFC3d SlipModel, is using friction angle. CoulombCriterion
if (shearForce.squaredNorm() > maxFs) {
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.squaredNorm() > maxFs) {
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*/.dot(shearForce) /*active force*/;
}
}

Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
2021-02-10
2021-02-12
 Jérôme Duriez (jduriez) said on 2021-02-12: #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/irreversible component of the increment in shear displacement.
(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 ;-) )