# en averaging for Law2_ScGeom_ViscElPhys_Basic, bug or incorrect documentation?

Asked by Daniel Kiracofe on 2020-05-28

The damping constant is computed at each contact in order to fulfill the normal restitution coefficient en=(en1 en2)/(en1+en2)

However, this does not appear to be what is actually implemented. It looks to me like the code actually does en = (en1+en2)/2

For example, if I test two materials, one with en=0.9 and the other with en=0.1, the formula stated in the documentation would result in en = 0.9 * 0.1 / (0.9 + 0.1) = 0.09. Example code with dropping a ball onto a wall is below. The expected rebound height of the ball is square root of coefficient of restitution times initial height. i.e. expect a rebound height of 1% of the initial height. The actual rebound height that I get is 0.25, which implies the coefficient of restitution is 0.5. That is consistent with an arithmetic average (0.9+0.1)/2.

So, my question: should it be considered as the documentation is the intended behavior and this is an issue with the code? or the code is the intended behavior and this is an issue with the documentation? or neither one is right, and it should be something else?

The behavior in the documentation does not make physical sense to me. imagine two materials both with en=1. You would end up with (1*1)/(1+1) = 0.5. That doesn't make sense. en=1 implies elastic behavior, but that's not what would happen. So to me, I would suggest to keep current arithmetic averaging and change the documentation to match.

In reality, of course, coefficient of restitution is not really defined for an individual material, but only really for pair of materials (and even then only empirically determined). I looked around to see if I could find any articles of someone suggesting one averaging method or the other, but I could not find any.

###################################################################

import matplotlib.pyplot as pyplot
qt.View() #open the controlling and visualization interfaces

box_x = 0.05
box_y = 0.05
box_z = 0.05

particle_dia = 0.005

young2 = 200e9
rho = 8230

mn = Vector3(0, box_y, 0)
mx = Vector3(box_x,2*box_y, box_z) #corners of the initial packing
thick = 2*particle_dia # the thickness of the walls

bigmx = (mx[0], 3 * mx[1], mx[2])

restitution_ball = 0.9
restitution_wall = 0.1

O.materials.append(ViscElMat(young=young2,poisson=0.7,frictionAngle=0, en = restitution_ball , et=1, density=rho , label='ve_ball'))
O.materials.append(ViscElMat(young=young2,poisson=0.7,frictionAngle=0, en = restitution_wall , et=1, density=rho , label='ve_wall'))

ball = sphere( (mx[0]/2, 2*mx[1] , mx[2]/2 ) , particle_dia/2, material='ve_ball' )
ballIds = O.bodies.append(ball)

walls=utils.aabbWalls([mn,bigmx],thickness=thick,oversizeFactor=1.5,material='ve_wall')
wallIds=O.bodies.append(walls)

#turn on gravity and let it settle
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()] ),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()] ),
NewtonIntegrator(gravity=(0,-9.81,0),damping=0.0),
]
O.dt=.05*PWaveTimeStep()

plot.addData(time=O.time , pos = (O.bodies[ballIds].state.pos[1]-box_y)/0.15 )

plot.plots={'time':('pos')}
plot.plot()

#O.run(-1, True)

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Bruno Chareyre (bruno-chareyre) said on 2020-05-28: #1

Hi,
Thanks for detailed explanation. I will not be able to compile this info right now personaly, but I have one comment.
The change following your previous question [1] is now merge in main branch.

A question is, if you are going to review/use that code together with its documentation, maybe you could suggest improvements directly via 'git' merge requests.
You can work on a separate branch to accumumate incremental fixes, then we discuss merge with ready.

Bruno

 Revision history for this message Daniel Kiracofe (kiracodl) said on 2020-05-28: #2

Yes, I am happy to make some suggested improvements via git merge requests. I'm not sure if I will find any more changes to this module or not, but I can make these changes now and if anything else related comes up later I can put that in too.

 Revision history for this message Robert Caulk (rcaulk) said on 2020-05-29: #3

> It looks to me like the code actually does en = (en1+en2)/2
>For example, if I test two materials...

Your empirical method seems to be confirmed by the source code ;) [1]

>The behavior in the documentation does not make physical sense to me.
>I looked around to see if I could find any articles of someone suggesting one averaging method or the other, but I could not find any.

I think I agree with your rationale but as I mentioned before, I am not up to date on that literature. The source code does reference an article "Pournin2001"[2]. However, I skimmed that article and it only defines kn, cn, kt, ct, while it doesn't seem to explicitly discuss the en for a contact (or maybe you have a bit more time to read deeper into it?).

Another, option is to email directly the original author, Sergei Dorofeenko <email address hidden>, to confirm the suspected documentation mistake.

[2]L. Pournin, Th. M. Liebling, A. Mocellin (2001), Molecular-dynamics force models for better control of energy dissipation in numerical simulations of dense granular media. Phys. Rev. E (65), pages 011302. DOI 10.1103/PhysRevE.65.011302

 Revision history for this message Bruno Chareyre (bruno-chareyre) said on 2020-05-29: #4

I think the equation therein have been modified a lot since Serguei Dorofenko, it is just not reflected in file authors.
It is sometimes difficult to answer "why is it averaged this way"?
"how to average" can be simply solutionless, then we pick a formula...
If there is rational for another averaging that's good to hear.
B