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

I use the simple-

[Ip2_

[Law2_

and setting NewtonIntegrato

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.

The source (pkg/dem/

How reliable is this code in Yade? Has someone double-checked that it does what it should?

Cheers, Václav

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- Yade Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Václav Šmilauer

- Solved:
- 2013-09-20

- Last query:
- 2013-09-20

- Last reply:
- 2013-09-18

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.

Anton Gladky (gladky-anton) said : | #2 |

Hi Vaclav,

probably this one:

http://

Anton

2013/9/18 Václav Šmilauer <email address hidden>:

> New question #235934 on Yade:

> https:/

>

> I use the simple-

>

> [Ip2_FrictMat_

> [Law2_ScGeom_

>

> and setting NewtonIntegrato

>

> 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.

>

> The source (pkg/dem/

>

> 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:/

> Post to : <email address hidden>

> Unsubscribe : https:/

> More help : https:/

Václav Šmilauer (eudoxos) said : | #3 |

The http://

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://

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://

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-

Vector3r normTemp = phys->normalForce - phys->normalVis

// 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->

if (normTemp.

phys-

phys-

}

else{

}

else if (includeAdhesion && phys->adhesionF

// *** limit viscous component to the max adhesive force ***

if (normTemp.

Real normVisc = phys->normalVis

phys-

phys-

}

// *** apply viscous component - in the presence of adhesion ***

else {phys->normalForce -= phys->normalVis

}

if (calcEnergy) {normDampDissip += phys->normalVis

}

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

Klaus Thoeni (klaus.thoeni) said : | #6 |

> 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://

@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:/

>

> 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://

>

> --

> You received this question notification because you are a direct

> subscriber of the question.

>

Yuxuan Wen (wenyuxuan) said : | #8 |

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(

in line 58. See https:/

But in Tsuji's paper, it is Kso = 2*sqrt(

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

Thank you.

Regards,

Yuxuan