different forces for sphere <-> facet(s) contact

Hello,
I was convinced that we answered the question.. If a sphere falls on the
facets boundary (i.e. on two facets at the same time), two interactions are
created and double force is applied on the particle (from two facets
instead of just one..). So in this specific case, the contact behaves like
it would have double stiffness (athought not exactly, I don't know why..).
HTH
Jan

2013/5/30 Alexander Eulitz [Eugen] <email address hidden>

> New question #229901 on Yade:
> https://answers.launchpad.net/yade/+question/229901
>
> Hi,
> we talked about this some time before, when I reported about the problem
> about "loosing spheres" in a simulation (I uploaded a video)
> Here is an example (corresponding code below). A sphere is falling caused
> by gravitation on the bottom of facetBox.
> It becomes apparent that contact forces (i.e. force of the sphere in z
> direction) are dependent on the position the sphere hits the box.
> If the sphere hits the box on a place where there two facets adjoin
> contact force will be higher than on any other place.
> What is the reason for this?
>
>
> #!/usr/bin/python
> # coding: utf-8
>
> from yade import pack, plot, geom, utils, timing
>
> import os
> import sys
> import time
>
> utils.readParamsFromTable(nr='NonBatch', height=0.05, r_sp=0.011980,
> m_sp=0.00351, contactLaw='HertzMindlin', E_sp=190e9, P_sp=0.30,
> Rho_sp=0.41, E_fl=190e9, P_fl=0.3, Rho_fl=0.41, density_fl=7835, dmp=0.2)
>
> from yade.params.table import *
> fileid=nr+str(height)+contactLaw
> if utils.runningInBatch(): fileid=O.tags['description']
>
> ##globale Variablen
> Messfrequenz=200000 #in Hz
> plotVirtPeriod=1.0/(Messfrequenz) #noch nicht eingebaut.
>
> #Density of the material [kg/m³]
> volume_sp=r_sp**3*4/3*math.pi
> density_sp=m_sp/volume_sp
> print 'Kugel:',volume_sp,'m³', density_sp
> print 'Boden:',density_fl,'kg/m³'
>
> O.materials.append(FrictMat(young=E_fl, poisson=P_fl,
> frictionAngle=Rho_fl, density=density_fl, label='mat_ground'))
> O.materials.append(FrictMat(young=E_sp, poisson=P_sp,
> frictionAngle=Rho_sp, density=density_sp, label='mat_spheres'))
>
> sp_left=O.bodies.append(utils.sphere([-0.15,0,height+r_sp],r_sp,color=[1,1,0],
> material='mat_spheres'))
> sp_middle=O.bodies.append(utils.sphere([0,0,height+r_sp],r_sp,color=[1,1,0],
> material='mat_spheres'))
> sp_right=O.bodies.append(utils.sphere([0.15,0,height+r_sp],r_sp,color=[1,1,0],
> material='mat_spheres'))
> #ground=O.bodies.append(utils.sphere([0,0,-r_sp*100],r_sp*100, fixed=True,
> color=[0,1,0], material='mat_ground'))
> #pot=O.bodies.append(utils.geom.facetCylinder((0,0,height/2),0.054,height,orientation=Quaternion((0,
> 0, 0),0),segmentsNumber=20,wallMask=2, color=(0,1,0), fixed=True,
> material='mat_ground'))
> pot=O.bodies.append(utils.geom.facetBox((0,0,height/2),(0.3,0.3,height/2),orientation=Quaternion((0,
> 0, 0),0),wallMask=16, color=(0,1,0), fixed=True, material='mat_ground'))
>
> print 'Programmstart am ',time.strftime("%d.%m.%Y um %H:%M:%S Uhr")
> print "Batch-Parameter"
> print O.tags['params']
>
> # Engines entsprechend Kontaktgesetz auswählen (quick&dirty)
> if contactLaw == 'HertzMindlin':
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
> [Ip2_FrictMat_FrictMat_MindlinPhys()],
> [Law2_ScGeom_MindlinPhys_Mindlin()]
>
> ),
> NewtonIntegrator(damping=dmp, gravity=(0,0,-9.81)),
> PyRunner(virtPeriod=plotVirtPeriod,command="myAddPlotData()",
> label='plotter'),
> PyRunner(virtPeriod=0.05,command="printTime()")
> ]
> elif contactLaw == 'HertzLinearShear':
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
> [Ip2_FrictMat_FrictMat_MindlinPhys()],
> [Law2_ScGeom_MindlinPhys_HertzWithLinearShear()]
>
> ),
> NewtonIntegrator(damping=dmp, gravity=(0,0,-9.81)),
> PyRunner(virtPeriod=plotVirtPeriod,command="myAddPlotData()",
> label='plotter'),
> PyRunner(virtPeriod=0.05,command="printTime()")
> ]
> elif contactLaw == 'Cundall':
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
> [Ip2_FrictMat_FrictMat_FrictPhys()],
> [Law2_ScGeom_FrictPhys_CundallStrack()]
>
> ),
> NewtonIntegrator(damping=dmp, gravity=(0,0,-9.81)),
> PyRunner(virtPeriod=plotVirtPeriod,command="myAddPlotData()",
> label='plotter'),
> PyRunner(virtPeriod=0.05,command="printTime()")
> ]
>
> O.dt=1*utils.PWaveTimeStep()
> print 'Timestep:',O.dt
>
> def myAddPlotData():
> Fmess_l=utils.sumForces([0], (0,0,1))
> Fmess_m=utils.sumForces([1], (0,0,1))
> Fmess_r=utils.sumForces([2], (0,0,1))
> force=plot.addData(t=O.time, Fnormal_left=Fmess_l,
> Fnormal_middle=Fmess_m, Fnormal_right=Fmess_r)
>
> def printTime():
> print O.time
>
>
> plot.plots={'t':('Fnormal_left', 'Fnormal_middle', 'Fnormal_right')}
>
>
> O.saveTmp()
>
> plot.plot()
> utils.waitIfBatch()
>
>
>
> --
> 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
>