Reference for coeff of restitution in Hertz-Mindlin with non-linear damping

Hi Vaclav,
I don't know if someone else than Chiara can give a precise answer. And unfortunately I suspect she is no longer reading yade's lists. You may contact her directly.
If you find an answer to you question, please let us know.

Hi Vaclav,

probably this one:

http://dx.doi.org/10.1016/0032-5910(92)88030-L

Anton

2013/9/18 Václav Šmilauer <email address hidden>:
> New question #235934 on Yade:
> https://answers.launchpad.net/yade/+question/235934
>
> I use the simple-scene-plot.py script to test the coefficient of restition for viscous contact laws, only replacing the two lines in the InteractionLoop with
>
> [Ip2_FrictMat_FrictMat_MindlinPhys(en=.5)],
> [Law2_ScGeom_MindlinPhys_Mindlin()],
>
> and setting NewtonIntegrator(damping=0.,...).
>
> The scene is a sphere in free fall onto a box, then rebounding, falling down again etc, until all energy is dissipated. According to the documentation and source code, *en* is the normal coefficient of restitution, so I expect to see the energy quartered (1/4=.5^2) after each rebound. The sphere is originally with v=0 at z=2, entering in contact with the box when z=1.5. Therefore, after one rebound, it should fly up to about 1.625 etc, eventually approaching 1.5.
>
> In the simulation, the sphere peaks at about z=1.55 after the first rebound. Considering that the equilibrium position will be at about z=1.45 due to still weight of the sphere in gravity, I still obtain coefficient of restitution of about sqrt((1.55-1.45)/(2-1.45))=.42, not 0.5.
>
> The source (pkg/dem/HertzMindlin.cpp, lines 89 and 320) references some (Tsuji, 1992) paper for computing the viscous coefficient from the coefficient of restitution (cn=alpha*sqrt(mbar), with alpha=-sqrt(5/6.)*2*log(en)/sqrt(log(en)^2)+M_PI^2)*sqrt(2*E*sqrt(R) and mbar=(m1*m2)/(m1+m2)). I am not able to identify this paper, can I ask for a full reference (e-mailing fulltext is highly appreciated)?
>
> How reliable is this code in Yade? Has someone double-checked that it does what it should?
>
> Cheers, Václav
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp

The http://people.ds.cam.ac.uk/jae1001/CUS/research/pfizer/Antypov_Elliott_EPL_2011.pdf paper gives the same equation under eqs (10) and (17) (though written a bit differently). I suspect that this paper is the source for the equation in Yade, since they cite Tsuji 1992 specifically as *not* giving analytical relationship for alpha(epsilon).

They discuss something just before the end of the paper ("Model choice"): the contact can be terminated when either the overlap reaches zero, or when the force is zero. With viscosity, there may be (unphysical) attractive force while there is still overlap, which might account for the extra dissipated energy.

I tried to limit (elastic+viscous) force to zero (so that there is no attraction due to viscosity). The plot is at http://flux.doxos.eu/sphere-wall-hertz.pdf : I get the energy balance almost exactly (max kinetic energy before the impact 5200J, after the impact, 3900J is dissipated by the viscous force, which makes exactly (5200-3900)/5200=0.25=0.5^2. The sphere rebounds from 1m to 0.27m, the still position is at -0.01315, giving coefficient of restitution of about sqrt(0.27+0.01315)=0.53. The 0.5 vs. 0.53 could be just my imprecise reading from plots, I assume.

I also tested collision time, which is given analytically in the paper above (eq (23)) and it agrees with simulation at different impact velocities (max. 1% difference).

I use Woo code now, but as far as I see, it should be equivalent to the code in Yade, except for the noAttraction flag mentioned above -- http://bazaar.launchpad.net/~eudoxos/woo/trunk/view/head:/pkg/dem/Hertz.cpp .

Hey Guys,

I am still here! Sorry for not replying earlier but I have recently been very busy and the good news is that I have successfully completed my PhD :-) (finally!). I take this occasion to say a Big Thanks to all of you for your friendship and support during this long journey which I greatly enjoyed, with all the lessons and gratifications that have come with it.

Now here is my reply to Vaclav (thanks for pointing this matter out).
You have spotted the problem since the normal force can become negative when viscous damping is applied. This should obviously be avoided as it is unphysical and below is the code that I have replaced in my local version of Yade that I have missed to commit (sorry about that). When adhesion is included, the code takes also care of this by limiting the normal force with the maximum adhesion force.

This issue can potentially affect all contact laws so ideally an option like this should be available even for the linear contact law when combined with some form of viscous damping.

Can somebody (Anton?) please replace the code below with what is currently in trunk in HertzMindlin.cpp, lines 344 to 345? Many thanks for this.

Let me know if there is any further question.

Cheers,
Chiara

 /**************************************/
 /* VISCOUS DAMPING (Normal direction) */
 /**************************************/
  // normal force must be updated here before we apply the Mohr-Coulomb criterion
 if (useDamping){ // get normal viscous component
  phys->normalViscous = cn*incidentVn;
  Vector3r normTemp = phys->normalForce - phys->normalViscous; // temporary normal force
  // viscous force should not exceed the value of current normal force, i.e. no attraction force should be permitted if particles are non-adhesive
  // if particles are adhesive, then fixed the viscous force at maximum equal to the adhesion force
  // *** enforce normal force to zero if no adhesion is permitted ***
  if (phys->adhesionForce==0.0 || !includeAdhesion){
   if (normTemp.dot(scg->normal)<0.0){
    phys->normalForce = Vector3r::Zero();
    phys->normalViscous = phys->normalViscous + normTemp; // normal viscous force is such that the total applied force is null - it is necessary to compute energy correctly!
   }
   else{phys->normalForce -= phys->normalViscous;}
  }
  else if (includeAdhesion && phys->adhesionForce!=0.0){
   // *** limit viscous component to the max adhesive force ***
   if (normTemp.dot(scg->normal)<0.0 && (phys->normalViscous.norm() > phys->adhesionForce) ){
    Real normVisc = phys->normalViscous.norm(); Vector3r normViscVector = phys->normalViscous/normVisc;
    phys->normalViscous = phys->adhesionForce*normViscVector;
    phys->normalForce -= phys->normalViscous;
   }
   // *** apply viscous component - in the presence of adhesion ***
   else {phys->normalForce -= phys->normalViscous;}
  }
  if (calcEnergy) {normDampDissip += phys->normalViscous.dot(incidentVn*dt);} // calc dissipation of energy due to normal damping
 }

Sounds good. Thank you Vaclav and Chiara for clarifications.
Now the only thing is to merge Vaclav's explanations and Chiara's code into an update of HM...
For the other viscous laws, it could also be fixed indeed. We will see.

Yet, I wonder if people do not calibrate the restitution vs. viscosity relation taking into account the tensile effects sometimes.

Cheers.

Bruno

> Yet, I wonder if people do not calibrate the restitution vs. viscosity
> relation taking into account the tensile effects sometimes.

At the end it depends what relation between the restitution coefficient and the
viscous parameter is used. There is even a paper on this topic [1].

Klaus

[1] http://link.springer.com/article/10.1007/s10035-007-0065-z

@Klaus
Yes, you are absolutely right and you refer to an excellent reference
(worth reading it for those interested in conctact models). However,
sometimes an analytical relationship is not possible to derive, such as in
the HM model.

@Bruno
At present I am not able to commit that part of the code but if somebody
follows my instruction it should be quick to do it.

All,
In my experience is not so critical to avoid tensile stresses due to
viscous damping as I have always found them to be negligible but it is a
good practice to neglect them entirely if possible. It is more of an issue
if we are modelling cohesive particles, as there you expect the cohesion to
match with your physical input.

Chiara

On 23 September 2013 11:16, Klaus Thoeni <
<email address hidden>> wrote:

> Question #235934 on Yade changed:
> https://answers.launchpad.net/yade/+question/235934
>
> Klaus Thoeni posted a new comment:
> > Yet, I wonder if people do not calibrate the restitution vs. viscosity
> > relation taking into account the tensile effects sometimes.
>
> At the end it depends what relation between the restitution coefficient
> and the
> viscous parameter is used. There is even a paper on this topic [1].
>
> Klaus
>
> [1] http://link.springer.com/article/10.1007/s10035-007-0065-z
>
> --
> You received this question notification because you are a direct
> subscriber of the question.
>

Hello Guys,

I just checked the source file of Hertz-Mindlin model and found that the shear stiffness ks is a little different with that in Tsuji's paper (1992). In the source file it is expressed as:
    Real Kso = 2*sqrt(4*R)*G/(2-V); // coefficient for shear stiffness"
in line 58. See https://github.com/yade/trunk/blob/master/pkg/dem/HertzMindlin.cpp
But in Tsuji's paper, it is Kso = 2*sqrt(2*R)*G/(2-V), so a factor of sqrt(2) exists between these 2 expressions.

Does it matter or is it just a typo? Should we correct it?

Thank you.

Regards,
Yuxuan