Understanding why sometimes does not occurs contact between Facet/Spheres using the same script

Asked by Marcus Moravia

Dear all,

I am starting a script to reach a specified porosity using REFD method. But right at the beginning, I have found some code behaviors that I still do not understand. In the script below:

1) If I use num=1 in the makeCloud there are no contacts between Facet/Sphere;
2) In the same script, if I activate the gravity=(0,0,-9,81) in NewtonIntegrator, there are contacts between Facet/Sphere;
3) Even without using gravity, but if I use, for example, 100 spheres (num=100 in the makeCloud), there are contacts between Facet/Sphere.

I'm a little unsure of these behaviors, could you tell me what's going on?

Thank you very much.

Marcus.

#############
## Script used ##
#############

from yade import utils,pack,geom,qt,plot

radius=.05
damp=0.2

O.materials.append(FrictMat(young=1e10,poisson=0.3,density=1e3,frictionAngle=radians(20),label='sphereMat'))
O.materials.append(FrictMat(young=1e10,poisson=0.3,density=1e3,frictionAngle=0,label='facetMat'))

O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=63,material='facetMat'))

sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=radius,num=100)
sp.toSimulation(material='sphereMat',color=(0,0,1))

triax=TriaxialStressController(
 maxMultiplier=0.01,
 finalMaxMultiplier=0.001,
 stressMask = 7,
 internalCompaction=True,
)

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(gravity=(0,0,0),damping=damp),
 PyRunner(command='radiusExpansion()',iterPeriod=200),
 PyRunner(command='contactForces()',iterPeriod=200),
]

O.dt=.5*PWaveTimeStep()

O.trackEnergy=True

def radiusExpansion():
 for n in range(0,len(O.bodies)):
  if isinstance(O.bodies[n].shape,Sphere):
   if O.bodies[n].shape.radius<0.6:
    O.bodies[n].shape.radius=O.bodies[n].shape.radius+0.0005

def contactForces():
 for i in O.interactions:
  if not i.isReal: continue
  cp = i.geom.contactPoint
  nf = i.phys.normalForce
  sf = i.phys.shearForce
  print cp

def fabricTensor():
 utils.plotDirections()

rr=yade.qt.Renderer()

plot.plots={'i':('unbalanced'),'i ':('sxx','syy','szz'),' i':('exx','eyy','ezz'),
 ' i ':(O.energy.keys,None,'Etot'),
}

yade.qt.View()
O.saveTmp()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:
Revision history for this message
Best Jan Stránský (honzik) said :
#1

Hi Marcus,
if you open the 3d view, switch off "shape" and switch on "bound", you will see that for num=1 you update radius of the sphere, but the bound is not updated. You can make the bound updated e.g. by setting
InsertionSortCollider(...,verletDist=0),
cheers
Jan

Revision history for this message
Marcus Moravia (mgmoravia) said :
#2

Thanks Jan Stránský, that solved my question.

Revision history for this message
Marcus Moravia (mgmoravia) said :
#3

Hi Jan,
Thank you very much for the trick of the switch off buttons. I could see exactly what was happening.
Cheers,
Marcus.