New question #691021 on Yade:
https://answers.launchpad.net/yade/+question/691021

The documentation for Law2_ScGeom_ViscElPhys_Basic 
(https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Law2_ScGeom_ViscElPhys_Basic)
 states the following:

  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
from yade import qt, plot
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),
        PyRunner(command='addPlotData()', iterPeriod=1000),
]
O.dt=.05*PWaveTimeStep()

def addPlotData():
        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) 




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

Reply via email to