float particles

Asked by littleFish

Dear all,
I have run a very simple script to prepare a packing with 10,000 particles. But when this process is finished with unbalancedForce lower than 0.01, I find there are more than 2,000 particles out of the total 10,000 particles without any interactions with other particles or walls ('float particles' for which O.bodies[i.id].intrs()==[]). I wonder if this is normal or there is something wrong with my simulation. The codes are attached below and thanks for your consideration.

from yade import pack, qt, plot
import time
import numpy as np

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

ts = time.time()
pressure = -1e5
r = 0.005
n = 10000
size0 = 0.3
mn,mx=Vector3(0,0,0),Vector3(size0,size0,size0)

## create materials for spheres and plates
O.materials.append(CohFrictMat(density=2650,young=6e8,poisson=.8,frictionAngle=0.5,isCohesive=True,normalCohesion=5e10,shearCohesion=5e10,momentRotationLaw=False,label='spheres'))
O.materials.append(CohFrictMat(density=0,young=8.8e10,poisson=.8,frictionAngle=0.,isCohesive=False,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()
sp.makeCloud(mn,mx,rMean=r,rRelFuzz=0.4,periodic=False,num=n,seed=1)
sp.toSimulation(material='spheres')

print(len(O.bodies))
############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
 ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
 ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
 #maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 #finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
 ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
 ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
 ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
 ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
 stressMask = 7,
 internalCompaction=False, # If true the confining pressure is generated by growing particles
        goal1 = pressure,
        goal2 = pressure,
        goal3 = pressure,
        max_vel = 1,
)

newton=NewtonIntegrator(gravity=(0,0,0), damping=.1)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),],
  [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionOnNewContacts=False,label='cohesiveIp'),],
  [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),]
 ),
 ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 #TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton,
        PyRunner(command='addPlotData()',iterPeriod=1000),
        PyRunner(command='check()',iterPeriod=1000),
        PyRunner(command='savePack()',iterPeriod=1000),
]

def addPlotData():

   pp = utils.porosity() #this is the porosity of the cell.
   ee = pp / (1-pp) #this is the void ratio of the 3D cell.
   a = [i for i in O.bodies if i.intrs()==[]]
   plot.addData(unbalanced=unbalancedForce(),i=O.iter,
                s11 = -getStress()[0,0],
                s22 = -getStress()[1,1],
                s33 = -getStress()[2,2],
                void = ee,
                Num = len(a)
   )

plot.plots={'i':('unbalanced',),'i ':('s11','s22','s33'),
   ' i ':('void',),' i':('Num')
}
# show the plot
plot.plot()

yade.qt.Controller(), yade.qt.View()

def savePack():
   O.save('packing_100kPa_Tmp.yade.gz')

def check():
  unb=unbalancedForce()
  if unb<0.01 and abs(pressure-triax.meanStress)/-pressure<0.001:
    print('Packing generated!')
    print('time: '+str((time.time()-ts)/60.)+'min')
    O.save('packing_100kPa.yade.gz')
    O.pause()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Bruno Chareyre
Solved:
Last query:
Last reply:
Revision history for this message
Robert Caulk (rcaulk) said :
#1

Thanks for posting a working script that demonstrates your question :-)

IINM you will always have some floaters when you use a uniform distribution of spheres like this, especially with a range associated with rRelFuzz=0.4 (very high!). 20% is very high. So you can decrease your floaters by 3 methods 1/ reducing the range (reducing rRelFuzz), 2/ reaching a specified porosity after you reach your desired confining pressure. Something like this (from examples/triax-tutorial[1]):

while triax.porosity>targetPorosity:
 # we decrease friction value and apply it to all the bodies and contacts
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 # while we run steps, triax will tend to grow particles as the packing
 # keeps shrinking as a consequence of decreasing friction. Consequently
 # porosity will decrease
 O.run(500,1)

3/ you can create "distant" by increasing the interactionRadius [2]. This will also decrease the number of floaters by creating a denser network. Warning: this will increase the number of interactions for networked particles too.

Cheers,

Robert

[1]https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py
[2]https://yade-dem.org/doc/user.html#creating-interactions

Revision history for this message
littleFish (littlefish) said :
#2

Thanks for your careful reply Robert!
I tried to set rRelFuzz=0 and the final packing still has about 1,000 floating particles out of 10,000. Since I need a fixed friction angle, it seems the only way is to increase aabbEnlargeFactor. Thanks again for your help though it's really strange that a simple packing has such a high proportion of floating particles.

Revision history for this message
Robert Caulk (rcaulk) said :
#3

>>it's really strange that a simple packing has such high proportion of floating particles.

How many floating particles *should* a randomly-placed single-sphere sized packing contain?

>>Since I need a fixed friction angle,

The process of finding a specified porosity ends with setting a fixed friction angle of your choice.

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

>it's really strange that a simple packing has such a high proportion of floating particles

It is actually not strange. Only organized packings leads to contacts on every particles. Others don't, and it is definitely not an issue. It has no relation with fixing a friction angle.

Bruno

Revision history for this message
littleFish (littlefish) said :
#5

Thanks Bruno Chareyre, that solved my question.