Yade

Validation of the Hertz-Mindlin model

Asked by Bettina Suhr

Dear all,

I found a nice publication for the verification of the Hertz-Mindlin contact model:
Chung, Y. & Ooi, J. Benchmark tests for verifying discrete element modelling codes at particle impact level Granular Matter, Springer-Verlag, 2011, 13, 643-656

There are 8 simple test defined for verification. For the first test (normal contact of two identical spheres) my result produced with Yade agrees with the analytical solution.
The next test is about normal contact between a sphere and a plane and here I get differences to the analytical solution.
I use the following script:

from yade import plot
import math
import sys
import numpy as np

O.materials.append(FrictMat(young=7.0e10, poisson=0.3, frictionAngle=np.arctan(0.0), density=2699.,label="alu"))
O.bodies.append([
      sphere(center=(0,0,0.11),radius=0.1, material='alu')
])
O.bodies.append(utils.wall(0, axis=2, sense=1))

#asign initial velocity
O.bodies[0].state.vel=(0,0,-0.2)

O.engines=[
  ForceResetter(),
  #use aabb for approx. contact detection
  InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(), Bo1_Facet_Aabb()]),
  InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_MindlinPhys()],#ELASTIC = No Damping
  [Law2_ScGeom_MindlinPhys_Mindlin()]
  ),
  NewtonIntegrator(damping=0.0, gravity=(0,0,0)),#NO DAMPING IN INTEGRATOR
  PyRunner(command='myMonitor()',iterPeriod=1)
]
O.dt=0.1*utils.PWaveTimeStep() # extra save value
print 'time step used: ', O.dt, ' critical time step size ', utils.PWaveTimeStep()
O.trackEnergy=True

t_start_contact=0

def myMonitor():
  global t_start_contact
  if len(O.interactions)>0 and O.interactions[0,1].geom:
    if np.abs(t_start_contact)<1e-10:
      t_start_contact=O.time
    for i in O.interactions:
       plot.addData(d=i.geom.penetrationDepth,Fn=np.abs(i.phys.normalForce[2]), t=O.time-t_start_contact)

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

# show the plot on the screen, and update while the simulation runs
plot.plot()

Do I do something wrong? Is this an error in Yade or is this a difference by intention?
Thanks for your help,
Bettina

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Bettina Suhr
Solved:
Last query:
Last reply:
Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#1

Interesting comparison!

For sphere vs. wall the wall is calculated as if it was a symmetric sphere. Not bad in general but in your case I guess it explains the difference, half-plane is more stiff than sphere.
The wal-sphere probably gives the same as sphere-sphere, no?

Honnestly most devs have been developping yade with bulk materials in mind, not really problems of one sphere impacting a plane. But it should be possible to fix that without a lot of difficulty.

Revision history for this message
Anton Gladky (gladky-anton) said : 		#2

It would also good to add those verified
tests to our "check"-scripts.

Bettina, it could be your good contribution into
Yade.

Anton

2015-07-28 17:46 GMT+02:00 Bruno Chareyre
<email address hidden>:
> Question #269739 on Yade changed:
> https://answers.launchpad.net/yade/+question/269739
>
> Status: Open => Answered
>
> Bruno Chareyre proposed the following answer:
> Interesting comparison!
>
> For sphere vs. wall the wall is calculated as if it was a symmetric sphere. Not bad in general but in your case I guess it explains the difference, half-plane is more stiff than sphere.
> The wal-sphere probably gives the same as sphere-sphere, no?
>
> Honnestly most devs have been developping yade with bulk materials in
> mind, not really problems of one sphere impacting a plane. But it should
> be possible to fix that without a lot of difficulty.
>
> --
> 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

Revision history for this message
Bettina Suhr (bettina-suhr) said : 		#3

Hi Bruno and Anton,

thanks for your comments.
@Bruno: You are right: the contact stiffness of the wall-sphere contact is equal to the case of symmetric spheres.
It is too low compared to the analytic solution. Can you explain me, what I need to change?
I was looking at Ig2_Wall_Sphere_ScGeom() and Ip2_FrictMat_FrictMat_MindlinPhys(), but I can't see any problem here. Or are there other functions which need to be changed?

Thanks,
Bettina

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

Speaking of wall-sphere, the definition of symmetric sphere radius is in [1].
This functor could optionaly define the radius of the virtual sphere in a different way.
Something to check is if analyticaly the 2-sphere equation goes to sphere-wall equation asymptotically when the radius of one sphere is increased (it should, no?). In such case you would only have to assign a very large value to the virtual radius (1e8*sphereRadius) and you would get the correct result without changing anything else.

Note that Box/Facet are more maintained than Walls, so maybe you could switch to them. A box with null thickness is a wall...

Bruno

[1] https://github.com/yade/trunk/blob/master/pkg/dem/Ig2_Facet_Sphere_ScGeom.cpp#L196

Revision history for this message
Bettina Suhr (bettina-suhr) said : 		#5

Dear Bruno,

thank you for you help. You are right, starting from the analytical solution for the contact of two spheres,
the analytical solution for sphere-wall contact is obtained when one radius goes to infinity.
What I found in the Yade source code is:

sphere-facet (Ig2_Facet_Sphere_ScGeom):
scm->radius1 = 2*sphereRadius;
scm->radius2 = sphereRadius;

sphere-wall contact (Ig2_Wall_Sphere_ScGeom):
ws->radius1=ws->radius2=radius; // do the same as for facet-sphere: wall's "radius" is the same as the sphere's radius

I changed in both function the radius of the wall to be big (1e8*radius of the sphere).

After eliminating another error in my script (my wall had the same material parameters as the sphere, where the wall was supposed to be rigid),
now I obtain the correct result for this test case.

Is the choice of the radius of a facet/wall meaningful for other contact laws? Or should it be changed?

Thanks again,
Bettina

Revision history for this message
Chareyre (bruno-chareyre-9) said : 		#6

Unfortunately this is not really helpful in general. Having boundary-sphere
more stiff than sphere-sphere can only increase heterogeneity near the
boundaries, hence maximizing boundary effects that we usually want to
avoid. It is interesting
as an option, not better in general. A proper way to give the option
without changing the default would be great.
Bruno

On 4 August 2015 at 17:46, Bettina Suhr <
<email address hidden>> wrote:

> Question #269739 on Yade changed:
> https://answers.launchpad.net/yade/+question/269739
>
> Status: Answered => Solved
>
> Bettina Suhr confirmed that the question is solved:
> Dear Bruno,
>
> thank you for you help. You are right, starting from the analytical
> solution for the contact of two spheres,
> the analytical solution for sphere-wall contact is obtained when one
> radius goes to infinity.
> What I found in the Yade source code is:
>
> sphere-facet (Ig2_Facet_Sphere_ScGeom):
> scm->radius1 = 2*sphereRadius;
> scm->radius2 = sphereRadius;
>
> sphere-wall contact (Ig2_Wall_Sphere_ScGeom):
> ws->radius1=ws->radius2=radius; // do the same as for facet-sphere: wall's
> "radius" is the same as the sphere's radius
>
> I changed in both function the radius of the wall to be big (1e8*radius
> of the sphere).
>
> After eliminating another error in my script (my wall had the same
> material parameters as the sphere, where the wall was supposed to be rigid),
> now I obtain the correct result for this test case.
>
> Is the choice of the radius of a facet/wall meaningful for other contact
> laws? Or should it be changed?
>
> Thanks again,
> Bettina
>
> --
> 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
>
>