Strange behavior of spheres colliding with a moving facet

def fill_box():
   """ fills the box with particles """
   def create_layer(N,r,h,offset=0):
      """ N -- number of particles in the layer
          r -- radius of particles
          h -- height of the layer
          offset -- initial offset """
      if offset==0:
         offset=2*r
      ni=nj=range(int(math.sqrt(N)))
      for i in ni:
         for j in nj:
            O.bodies.append(utils.sphere((offset+j*r*3,offset+i*r*3,h),r))#,material='particles'))

   for i in range(1):
      create_layer(16,0.008,0.01+i*0.018,offset=0.01+float(random.randint(10,20))/2500)

def shake_table():
   """ move the table according to the sine law (period T=0.5) """
   global pos_facet1
   global pos_facet2
   h=0.02*math.sin(2*pi*O.time)
   pos_facet1[2]=h
   pos_facet2[2]=h
   O.bodies[8].state.pos=pos_facet1
   O.bodies[9].state.pos=pos_facet2

# handling engines
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      # handle sphere+sphere and facet+sphere collisions
      [Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_L3Geom_FrictPhys_ElPerfPl()]
   ),
   GravityEngine(gravity=(0,0,-9.81)),
   NewtonIntegrator(damping=0.04),
   PyRunner(command='shake_table()',virtPeriod=O.dt),
]
# time step
O.dt=1e-06

# MAIN SCENARIO

# Create a box with no roof
O.bodies.append(utils.facetBox((0.05,0.05,0.06),(.05,.05,.1),wallMask=15))
O.materials[0].density=1500
# Create a table
O.bodies.append(utils.facet([(0,0,0),(0,0.1,0),(0.1,0.1,0)],color=(1,0,0)))
O.bodies.append(utils.facet([(0,0,0),(0.1,0,0),(0.1,0.1,0)],color=(1,0,0)))
# copy their coords into variables
pos_facet1=O.bodies[8].state.pos
pos_facet2=O.bodies[9].state.pos
# fill the box with spheres
fill_box()

O.saveTmp()

Try to decrease time step (O.dt = 1e-8 or 1e-7) when facet comes in
contact with balls.

Christian

Zitat von vadim <email address hidden>:

> New question #171223 on Yade:
> https://answers.launchpad.net/yade/+question/171223
>
> Hello to everyone! I use Yade 1+2922+27~natty1
> The model consists of spheres, a moving table (two facets) and walls
> around it. Initially spheres (almost) lie on the table. When the
> table begins moving some spheres start moving with incredibly hight
> velocities and act very strange. It seems that collision
> sphere/facet is not handled properly. Am I doing something wrong?
> Thank you
> Vadim
>
> --
> 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
>

Since the table is shaking, I guess relatively large velocities can be
expected.
So... I turned off shaking, and then, indeed I see a bug.
Actually, L3-based functors are experimental and no longer maintained.
You would get more help, and probably less bugs, by using ScGeom functors.

1. Time step = 1e-6:

Yade [18]: O.run(20000)
Yade [19]: O.time
      -> [19]: 0.019999999999999348
Yade [20]: O.bodies[15].state.vel
      -> [20]: Vector3(0,0,3.0782119911613148)

2. Time step = 1e-7:

Yade [12]:O.run(200000)
Yade [13]: O.time
      -> [13]: 0.019999999999935486
Yade [14]: O.bodies[15].state.vel
      -> [14]: Vector3(0,0,3.0782980240123341)

3. Time step = 1e-8:

Yade [15]: O.run(2000000)
Yade [16]: O.time
      -> [16]: 0.020000000000648145
Yade [17]: O.bodies[15].state.vel
      -> [17]: Vector3(0,0,3.0783555450742548)

3 meters per second and time step doesn't seem to make changes! How can it be if the table is moving with amplitude of 2 mm and period of 1 sec? (0.02*math.sin(2*pi*O.time))

The bug is the same with ScGeom. It's not linked to the geometry type
(or let's say, L3 duplicated the bug).
You probably pointed out a very interesting bug, that is maybe related
to other bugs in facet-spheres interaction (e.g.
https://bugs.launchpad.net/yade/+bug/585898). Thanks for that.
If you want to investigate more, please report anything interesting
about the bug here: https://bugs.launchpad.net/yade/+bug/850864
If you want correct results, you can use box instead of facets.

I'm not expert in how impacts and vibrations are combined, but I can see
something is wrong even without the vibrations, so... no point vibrating
the bug ;)

The sphere is touching both facets.
Perhaps, it gets 2 forces of the same value instead of one.

Anton

I changed the model and made a "table" from just one facet:

===================================================
def fill_box():
   """ fills the box with particles """
   def create_layer(N,r,h,offset=0):
      """ N -- number of particles in the layer
          r -- radius of particles
          h -- height of the layer
          offset -- initial offset """
      if offset==0:
         offset=2*r
      ni=nj=range(int(math.sqrt(N)))
      for i in ni:
         for j in nj:
            O.bodies.append(utils.sphere((offset+j*r*3,offset+i*r*3,h),r))

   for i in range(1):
      create_layer(16,0.008,0.01+i*0.018,offset=0.01+float(random.randint(10,20))/2500)

def shake_table():
   """ move the table according to the sine law (period T=0.5) """
   global pos_table
   h=0.02*math.sin(2*pi*O.time*2)
   pos_table[2]=h
   O.bodies[8].state.pos=pos_table

# handling engines
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      # handle sphere+sphere and facet+sphere collisions
      [Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_L3Geom_FrictPhys_ElPerfPl()]
   ),
   GravityEngine(gravity=(0,0,-9.81)),
   NewtonIntegrator(damping=0.04),
   PyRunner(command='shake_table()',virtPeriod=O.dt),
]
# time step
O.dt=1e-08

# MAIN SCENARIO

# Create a box with no roof
O.bodies.append(utils.facetBox((0.05,0.05,0.06),(.05,.05,.1),wallMask=15))
O.materials[0].density=1500
# Create a table
O.bodies.append(utils.facet([(0.15,-0.1,0),(0.15,0.2,0),(-0.15,0.05,0)],color=(1,0,0)))
# copy coordinates of the table into a global variable
pos_table=O.bodies[8].state.pos
# fill the box with spheres
fill_box()

O.saveTmp()
===================================================

Let's compare new results with the previous:

1. Old, two facet table, time step dt=1e-8:

Yade [15]: O.run(2000000)
Yade [16]: O.time
       -> [16]: 0.020000000000648145
Yade [17]: O.bodies[15].state.vel
       -> [17]: Vector3(0,0,3.0783555450742548)

2. New, one facet table, time step dt=1e-8:

Yade [1]: O.run(2000000)
Yade [2]: O.time
      -> [2]: 0.020000000000648145
Yade [3]: O.bodies[15].state.vel
      -> [3]: Vector3(0,0,3.6318959155341539)

Conclusions. Number of facets doesn't matter. Yade doesn't work out sphere/facet interactions properly. Maybe I can use a sphere with huge radius instead of the facet?

PS. At this time all spheres jump with the same speed

> Conclusions. Number of facets doesn't matter. Yade doesn't work out
> sphere/facet interactions properly. Maybe I can use a sphere with huge
> radius instead of the facet?

As suggested before, you could use boxes, but actually your conclusion
does not hold, the problem is in your script.
Bodies motion should always be defined in terms of prescribed velocity
(or using motion Engines for the most complicated cases).
Here, you are changing the position of the shaking table, which breaks
the current implementation of verlet distance in the collider. The
result wuld be the same with facet, box, or a big sphere.

Setting verletDist=0 fixes the problem, but still it is not good:
viscous laws for instance would suffer from this pos= in the script.
Prescibe velocity, then it will work like a charm, and you can keep
using verlet>0 (default), which will speed up the simulations.

Thank you, Chareyre! It helped

But the behavior is still strange for me. I remade model a little bit: there is a moving table and just one sphere. The sphere falls and hits the table at step 362 150. Interaction 'lasts' until 362 950. I was going to find out how speed of the sphere is calculated...

------------------------------------------------------------------------------
Yade [167]: O.run(362150)
Yade [168]: pvel()

sphere.vel = Vector3(0,0,-1.6687395853455125)
table.vel = Vector3(0,0,0.80683491142746167)

Yade [169]: O.run(362950-362150)
Yade [170]: pvel()

sphere.vel = Vector3(0,0,3.1912989954315774)
table.vel = Vector3(0,0,0.79939792265750198)
------------------------------------------------------------------------------

My question is: if speed of the sphere before interaction is 1.67 and speed of the table is 0.81. And 1.67+0.81 = 2.48. So how can speed of the sphere be 3.19 after interaction? It should be less then 2.48 because of loss but it's greater! Probably I don't understand something

Vadim

Vadim,

I think the result is physically correct. If you think in a moving frame of reference, you get incident velocity = 2.48 and post-impact velocity 3.19-0.81=2.38.
I guess you have a little bit of damping that can explain the speed decrease?

Chareyre,

If the frame of reference is moving (with table, right?), the "post-impact" velocity has to be 2.48, not 3.19. But what we have here?

1. Before interaction
     sphere velocity: -1.67
     table velocity: 0.81

2. After interaction
     sphere velocity: 3.19 (why not 0.81 - (-1.67) = 2.48 ?)
     table velocity: 0.81

The sphere got some extra energy

Yes, frame moving with table. Intuitively, I would think you have to define post-impact velocity in this moving frame as well.
3.19 is the velocity in absolute coordinates. In the moving frame it gives 3.19-0.81=2.38, no?
Extra energy is expected I guess, since the table is impacting the sphere.

Hello,

I think there is a time step problem. If the time step is to big, the overlap between the sphere and the facet becomes too big within one step. So it is possible, that the sphere gets a very big contact force and bounces back much faster, than it hits. I created a very small script for you, where you can see the effect:

#!/usr/bin/python
# -*- coding: utf-8 -*-

from yade import qt

#define material properties:
rho = 1000 #density
fric = 2 #friction coefficient
fric_angle = atan(fric) #friction angle in radian
stiff_n = 1e8
stiff_s = 1e8
c_n = 0.0 #normal viscous damping
c_s = 0.0 #shear viscous damping
l_d = 0.7 #local damping

#define a material:
MyMaterial = O.materials.append(ViscElMat(kn=stiff_n, ks=stiff_s, \
cn=c_n, cs=c_s, density=rho, frictionAngle=fric_angle))

#create a sphere:
O.bodies.append(utils.sphere(center=(0,0,1), radius=0.5, material=MyMaterial))

#create a plane facet:
O.bodies.append(geom.facetBox(center=(0,0,0), extents=(1,1,0), fixed=True, material=MyMaterial))

#define engines:
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]), #Axis-aligned bounding box
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom()], #Interaction geometry
  [Ip2_ViscElMat_ViscElMat_ViscElPhys()], #Interaction physics
  [Law2_ScGeom_ViscElPhys_Basic()], #Contact law
 ),
 GravityEngine(gravity=(0,0,-9.81)),
 NewtonIntegrator(damping=l_d),
 #qt.SnapshotEngine(fileBase='/home/me/YADE/my-yade-projects/sphere-and-facet/snaps/test-',label='snapper',iterPeriod=10),
]

#show geometry:
v=qt.View(); v.axes=True
v.viewDir=Vector3(0,1,-.05)
v.eyePosition=Vector3(0,-5,.5)

O.dt = 1e-3
O.run(1000)

#O.dt = 1e-6
#O.run(1000000)

#raw_input('press ENTER')
#utils.makeVideo(snapper.snapshots,out="video-sphere-and-facet.avi",fps=24)

#END of script##############################################

If you first run it with O.dt = 1e-3, the sphere bounces back very fast.
If you run it with O.dt = 1e-6, it seems to behave normal (bzr2877).
The same behavior should occur with using a fixed sphere instead of a facet (I did not test it)!

@Bruno:
I detected an bug in your new collider version!
If I run the script with O.dt = 1e-6 with bzr2877 the collision is ok, but with bzr2921 the sphere bounces back very fast, like in the O.dt = 1e-3 case.

I hope it helps,

Christian

Sorry,

My last comment is not true, I forgot to save python script before running it.

So there is no bug! New collider works well!

Christian

Yes, clearly, if dt is too large the impact will not be described correctly and it can create energy. This is a feature of explicit DEM, not a bug. You have to consider incident velocities when defining dt, independently of the numerical stability condition.
However, you will get the 3.19 velocity even with very small timesteps I think (and it is correct, do you agree?).

Thanks for testing the branch! :)

>1. Before interaction
> sphere velocity: -1.67
> table velocity: 0.81
>
>2. After interaction
> sphere velocity: 3.19 (why not 0.81 - (-1.67) = 2.48 ?)
> table velocity: 0.81

I think the velocity after interaction in absolute coordinates should be 2.48 (=1.67 + 0.81). With damping=0.04 it should be even a bit lower.

My idea was that the difference between 2.48 and 3.19 could occur, when the time step is too big.
I dont know if vadim tested this ...

Christian

Chareyre,
You are making the issue more complicated. Let's just think in absolute frame of reference

Chrisian,
I agree with you. The time step I used: 1e-6

Conclusions. I'll use YADE for my thesis. Even if there are some bugs, no one knows it because everyone uses MATLAB ;)

Vadim,
I'm ok to think in any frame of reference, but for now I don't see any evidence that something is wrong.
In absolute frame there is no reason to sum 1.67 + 0.81. The velocity is -1.67, then 3.19, so what?

The reasoning I suggest in moving frame is the only rational way I can find to predict a post-impact velocity, and it's confirmed by Yade, which is very unlikely to be wrong on such simple mass-spring problem.

Yes, of course. Now I get it.
You are right, everything is physically correct.
I absolutely agree with you.

Thanks Chareyre, that solved my question.