New question #695457 on Yade:
https://answers.launchpad.net/yade/+question/695457
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.
Thanks in advance,
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*/;
}
}
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
