Strange behaviour of Ek varying restitution coefficients

Asked by Alessandro Finazzi on 2019-04-09

Hi everybody,

I have to simulate the behaviour of some spheres in a moving cylinder varying the tangential and normal restitution coefficients of the material. What seems strange to me is that I get a kinetic Energy increasing with lower restitution coefficients.
From the definition of coefficient of restitution I would expect the contrary.

I am wondering if I'm making a silly conceptual mistake or if there is an error in my code.
Do you see any error in the following code?
I work with yadedaily, xenial version.
Many thanks

from yade import pack,utils,plot,ymport,export,qt,timing

## Physical Parameters
density=7950
cyl_r = 12.7e-3
cyl_h = 24.6e-3
sp_p_r = 1.5e-3
young = 200e9
poisson = 0.305

## Parameters to be varied
#en = 0.1
#es = 0.1
en = 0.78
es = 0.78

## Creation Cylinder
facetMat=O.materials.append(ViscElMat(density=density,young=young, poisson=poisson,en=en,et=es))
Cylinder=O.bodies.append(geom.facetCylinder((0,0,cyl_h/2),cyl_r,height=cyl_h,segmentsNumber=15,wallMask=7,material=facetMat))

## Creation Sphere
sphereMat=O.materials.append(ViscElMat(density=density,young=young,poisson=poisson,en=en,et=es))

s=utils.sphere((0,0,12e-3),radius=3e-3,material=sphereMat)
O.bodies.append(s)

sp=pack.SpherePack()
sp.makeCloud((-9e-3,-9e-3,1.5e-3),(9e-3,9e-3,8e-3),rMean=sp_p_r,num=20)
sp.toSimulation(material=sphereMat)

## Engine
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
      [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
      [Law2_ScGeom_ViscElPhys_Basic()]
   ),
   GlobalStiffnessTimeStepper(timeStepUpdateInterval=1,label='tss',timestepSafetyCoefficient=0.8,viscEl=True),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.2),
   PyRunner(command='checkUnbalanced()',realPeriod=1,label='checker')
   ]

def checkUnbalanced():
 if O.iter<100000: return
        O.engines=O.engines+[HarmonicMotionEngine(A=[0,0,6e-3],f=[0,0,60.0],fi=[0,0,pi],ids=Cylinder), PyRunner(command='addPlotData()',realPeriod=1)]
 checker.command='nothing()'

def nothing():
        return

def addPlotData():
 Ek=kineticEnergy()
 plot.addData(i=O.iter,t=O.time,Ninter=O.interactions.countReal(),total=O.energy.total(),Ek=kineticEnergy(),**O.energy)

O.trackEnergy=True
O.saveTmp()

plot.plots={'t':('Ek',)}

plot.plot()

qt.Controller()
qt.View()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Chareyre
Solved:
2019-04-09
Last query:
2019-04-09
Last reply:
2019-04-09
Best Chareyre (bruno-chareyre-9) said : #1

Hi, first of all you could have inconsistent restitution due to collisional
error (late detection of contact). Try decreasing the timestep and see if
it's better.

Second, contact viscosity is not very accurately integrated with centered
scheme, again decreasing timestep shoul help.

Finally, Cundall damping should be null for consistent results.

Le mar. 9 avr. 2019 18:08, Alessandro Finazzi <
<email address hidden>> a écrit :

> New question #680145 on Yade:
> https://answers.launchpad.net/yade/+question/680145
>
> Hi everybody,
>
> I have to simulate the behaviour of some spheres in a moving cylinder
> varying the tangential and normal restitution coefficients of the material.
> What seems strange to me is that I get a kinetic Energy increasing with
> lower restitution coefficients.
> >From the definition of coefficient of restitution I would expect the
> contrary.
>
> I am wondering if I'm making a silly conceptual mistake or if there is an
> error in my code.
> Do you see any error in the following code?
> I work with yadedaily, xenial version.
> Many thanks
>
>
>
> from yade import pack,utils,plot,ymport,export,qt,timing
>
> ## Physical Parameters
> density=7950
> cyl_r = 12.7e-3
> cyl_h = 24.6e-3
> sp_p_r = 1.5e-3
> young = 200e9
> poisson = 0.305
>
> ## Parameters to be varied
> #en = 0.1
> #es = 0.1
> en = 0.78
> es = 0.78
>
> ## Creation Cylinder
> facetMat=O.materials.append(ViscElMat(density=density,young=young,
> poisson=poisson,en=en,et=es))
>
> Cylinder=O.bodies.append(geom.facetCylinder((0,0,cyl_h/2),cyl_r,height=cyl_h,segmentsNumber=15,wallMask=7,material=facetMat))
>
> ## Creation Sphere
>
> sphereMat=O.materials.append(ViscElMat(density=density,young=young,poisson=poisson,en=en,et=es))
>
> s=utils.sphere((0,0,12e-3),radius=3e-3,material=sphereMat)
> O.bodies.append(s)
>
> sp=pack.SpherePack()
> sp.makeCloud((-9e-3,-9e-3,1.5e-3),(9e-3,9e-3,8e-3),rMean=sp_p_r,num=20)
> sp.toSimulation(material=sphereMat)
>
> ## Engine
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
> [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
> [Law2_ScGeom_ViscElPhys_Basic()]
> ),
>
> GlobalStiffnessTimeStepper(timeStepUpdateInterval=1,label='tss',timestepSafetyCoefficient=0.8,viscEl=True),
> NewtonIntegrator(gravity=(0,0,-9.81),damping=0.2),
> PyRunner(command='checkUnbalanced()',realPeriod=1,label='checker')
> ]
>
> def checkUnbalanced():
> if O.iter<100000: return
>
> O.engines=O.engines+[HarmonicMotionEngine(A=[0,0,6e-3],f=[0,0,60.0],fi=[0,0,pi],ids=Cylinder),
> PyRunner(command='addPlotData()',realPeriod=1)]
> checker.command='nothing()'
>
> def nothing():
> return
>
> def addPlotData():
> Ek=kineticEnergy()
>
> plot.addData(i=O.iter,t=O.time,Ninter=O.interactions.countReal(),total=O.energy.total(),Ek=kineticEnergy(),**O.energy)
>
>
> O.trackEnergy=True
> O.saveTmp()
>
> plot.plots={'t':('Ek',)}
>
> plot.plot()
>
> qt.Controller()
> qt.View()
>
> --
> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>
>
>

Hi Bruno,

it seems that, as you suggested, decreasing the time step is the right solution.
With a time step that is a tenth of my original one, I can get plausible results even with the smallest coefficients I tried (en=et=0.01).
Thank you very much, also for the suggestion about the damping.

Thanks Chareyre, that solved my question.