How to count number of particles in a fragment, connected by cohesive bonds

Asked by VG on 2016-04-29

Hello,

I am trying to model crushing of aggregates of particles and look at the resulting fragment size distribution. The model consists of an assembly of spheres (each sphere is composed of constituent smaller particles) between the two plates. The top plate is subjected to a constant force and a velocity in lateral direction to induce shear.

As the simulation progresses, the cohesive bonds between the different particles start breaking resulting into individual particles and fragments of particles. How can I measure the number of such individual particles, fragments and the number of particles in each fragment ?
I want to do this to obtain the particle size distribution of the crushed sample.

Here is the minimal working example script:

from yade import pack

#############################################################################
# Set up run
#############################################################################
run_name="PeriodicCohesive_MWE"
data_root_dir="."

#############################################################################
# Materials
#############################################################################
plate_material=FrictMat(
    young=200e9
   ,poisson=0.3
   ,density=8000.
   ,frictionAngle=radians(30)
   ,label='plate_mat')

O.materials.append(plate_material)

sample_material=CohFrictMat(
    young=4e9
   ,poisson=0.25
   ,density=1400
   ,frictionAngle=radians(30)
   ,normalCohesion=1e8*1.2
   ,shearCohesion=.4e8*1.2
   ,momentRotationLaw=True
   ,label='sample_mat')
O.materials.append(sample_material)

#############################################################################
# Component dimensions and operating condition
#############################################################################
# Granular material dimension
sample_diameter=2e-4
sample_radius=sample_diameter/2.0
# Sub-particle dimension
particle_diameter=74.e-6
particle_radius=particle_diameter/2.

#############################################################################
# grinding plate dimension
#############################################################################

rotvel=2./3.*pi*(1.5+0.5)*.254

#############################################################################
# Periodic Geometry
#############################################################################

# Set up periodic boundary conditions
O.periodic=True
xExt=4*sample_diameter
yExt=3.*sample_diameter*2 #to block the periodicity in y direction
zExt=xExt
xLim=xExt
yLim=yExt/4
zLim=zExt
O.cell.hSize=Matrix3(
  xExt, 0, 0,
  0, yExt, 0,
  0, 0, zExt)

length=xExt
height=yExt
width=zExt

# Top and bottom plate thickness
thickness=0.1*height

bottomBoxes = []
for ix in (0,1,2):
 for iz in (0,1,2):
   bottomBoxes.append(box( # create 3x3 boxes with 1/3 cell size
   center=(xExt/6.*(1+2*ix),yLim - thickness/2.0,zExt/6.*(1+2*iz))
   ,extents=(xExt/6.,thickness/2.0,zExt/6.)
   ,wire=False
   ,material='plate_mat'
   ))

bottom_id,bottom_ids = O.bodies.appendClumped(bottomBoxes) # bottom_id is the clump id,

O.bodies[bottom_id].state.blockedDOFs='xyzXYZ'

#############################################################################
# Particle Packing
#############################################################################

min_corner= (0,yLim,0)
max_corner= (xLim, yExt-yLim, zLim)

sp=pack.SpherePack()
sp.makeCloud( min_corner,max_corner, rMean=sample_radius, periodic=False)

print "Generated ",len(sp)," aggregates"

###########################################
# Sample
###########################################
for s in sp:
  sphere=pack.inSphere((s[0][0],s[0][1],s[0][2]),s[1])
  sp1=pack.randomDensePack(
      sphere
     ,spheresInCell=2000
     ,radius=particle_radius
     ,memoizeDb='/tmp/triaxPackCache.sqlite'
     ,returnSpherePack=True
     )

  sp1.toSimulation(material='sample_mat',color=(0.9,0.8,0.6))
  print 'Generated ',len(sp1),' particles'

Gl1_Sphere(stripes=True)

#########
 # Top
#########
topBoxes = []
for ix in (0,1,2):
 for iz in (0,1,2):
   topBoxes.append(box( # create 3x3 boxes with 1/3 cell size
   center=(xExt/6.*(1+2*ix),yExt - yLim + thickness/2.0,zExt/6.*(1+2*iz))
   ,extents=(xExt/6.,thickness/2.0,zExt/6.)
   ,wire=False
   ,material='plate_mat'
   ))

top_id,top_ids = O.bodies.appendClumped(topBoxes) # top_id is the clump id,

O.bodies[top_id].state.blockedDOFs='xzXYZ'

plate_downforce=-0.036

O.forces.addF(top_id,(0,plate_downforce,0),permanent=True)

O.bodies[top_id].state.vel[0]= rotvel

#############################################################################
# Run the simulation
#############################################################################

O.dt=0.5*PWaveTimeStep()

O.engines=[
  ForceResetter(),
  InsertionSortCollider([
    Bo1_Sphere_Aabb()
   ,Bo1_Box_Aabb()
   ], allowBiggerThanPeriod=True
   ),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom6D()
    ,Ig2_Box_Sphere_ScGeom6D()
    ],
    [Ip2_FrictMat_FrictMat_FrictPhys()
     ,Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,label="cohesiveIp")
    ],
    [Law2_ScGeom_FrictPhys_CundallStrack()
     ,Law2_ScGeom6D_CohFrictPhys_CohesionMoment()
    ]
  ), # End InteractionLoop
  NewtonIntegrator(damping=0.8,gravity=(0.,0.,0.)),

]

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2016-04-29
Last reply:
2016-05-02
Jan Stránský (honzik) said : #1

Hi Varun,
the solution should not be difficult. To test it, please modify your MWE to
be actually working, i.e. to contain correct interactions :-)
thanks
Jan

2016-04-30 2:02 GMT+02:00 VG <email address hidden>:

> New question #292911 on Yade:
> https://answers.launchpad.net/yade/+question/292911
>
> Hello,
>
> I am trying to model crushing of aggregates of particles and look at the
> resulting fragment size distribution. The model consists of an assembly of
> spheres (each sphere is composed of constituent smaller particles) between
> the two plates. The top plate is subjected to a constant force and a
> velocity in lateral direction to induce shear.
>
> As the simulation progresses, the cohesive bonds between the different
> particles start breaking resulting into individual particles and fragments
> of particles. How can I measure the number of such individual particles,
> fragments and the number of particles in each fragment ?
> I want to do this to obtain the particle size distribution of the crushed
> sample.
>
> Here is the minimal working example script:
>
>
> from yade import pack
>
>
> #############################################################################
> # Set up run
>
> #############################################################################
> run_name="PeriodicCohesive_MWE"
> data_root_dir="."
>
>
>
>
> #############################################################################
> # Materials
>
> #############################################################################
> plate_material=FrictMat(
> young=200e9
> ,poisson=0.3
> ,density=8000.
> ,frictionAngle=radians(30)
> ,label='plate_mat')
>
> O.materials.append(plate_material)
>
>
> sample_material=CohFrictMat(
> young=4e9
> ,poisson=0.25
> ,density=1400
> ,frictionAngle=radians(30)
> ,normalCohesion=1e8*1.2
> ,shearCohesion=.4e8*1.2
> ,momentRotationLaw=True
> ,label='sample_mat')
> O.materials.append(sample_material)
>
>
>
>
> #############################################################################
> # Component dimensions and operating condition
>
> #############################################################################
> # Granular material dimension
> sample_diameter=2e-4
> sample_radius=sample_diameter/2.0
> # Sub-particle dimension
> particle_diameter=74.e-6
> particle_radius=particle_diameter/2.
>
>
> #############################################################################
> # grinding plate dimension
>
> #############################################################################
>
> rotvel=2./3.*pi*(1.5+0.5)*.254
>
>
> #############################################################################
> # Periodic Geometry
>
> #############################################################################
>
>
> # Set up periodic boundary conditions
> O.periodic=True
> xExt=4*sample_diameter
> yExt=3.*sample_diameter*2 #to block the periodicity in y direction
> zExt=xExt
> xLim=xExt
> yLim=yExt/4
> zLim=zExt
> O.cell.hSize=Matrix3(
> xExt, 0, 0,
> 0, yExt, 0,
> 0, 0, zExt)
>
>
> length=xExt
> height=yExt
> width=zExt
>
> # Top and bottom plate thickness
> thickness=0.1*height
>
>
>
>
>
> bottomBoxes = []
> for ix in (0,1,2):
> for iz in (0,1,2):
> bottomBoxes.append(box( # create 3x3 boxes with 1/3 cell size
> center=(xExt/6.*(1+2*ix),yLim - thickness/2.0,zExt/6.*(1+2*iz))
> ,extents=(xExt/6.,thickness/2.0,zExt/6.)
> ,wire=False
> ,material='plate_mat'
> ))
>
> bottom_id,bottom_ids = O.bodies.appendClumped(bottomBoxes) # bottom_id is
> the clump id,
>
>
> O.bodies[bottom_id].state.blockedDOFs='xyzXYZ'
>
>
>
>
> #############################################################################
> # Particle Packing
>
> #############################################################################
>
> min_corner= (0,yLim,0)
> max_corner= (xLim, yExt-yLim, zLim)
>
> sp=pack.SpherePack()
> sp.makeCloud( min_corner,max_corner, rMean=sample_radius, periodic=False)
>
> print "Generated ",len(sp)," aggregates"
>
> ###########################################
> # Sample
> ###########################################
> for s in sp:
> sphere=pack.inSphere((s[0][0],s[0][1],s[0][2]),s[1])
> sp1=pack.randomDensePack(
> sphere
> ,spheresInCell=2000
> ,radius=particle_radius
> ,memoizeDb='/tmp/triaxPackCache.sqlite'
> ,returnSpherePack=True
> )
>
> sp1.toSimulation(material='sample_mat',color=(0.9,0.8,0.6))
> print 'Generated ',len(sp1),' particles'
>
> Gl1_Sphere(stripes=True)
>
> #########
> # Top
> #########
> topBoxes = []
> for ix in (0,1,2):
> for iz in (0,1,2):
> topBoxes.append(box( # create 3x3 boxes with 1/3 cell size
> center=(xExt/6.*(1+2*ix),yExt - yLim + thickness/2.0,zExt/6.*(1+2*iz))
> ,extents=(xExt/6.,thickness/2.0,zExt/6.)
> ,wire=False
> ,material='plate_mat'
> ))
>
> top_id,top_ids = O.bodies.appendClumped(topBoxes) # top_id is the clump id,
>
>
> O.bodies[top_id].state.blockedDOFs='xzXYZ'
>
>
>
> plate_downforce=-0.036
>
> O.forces.addF(top_id,(0,plate_downforce,0),permanent=True)
>
> O.bodies[top_id].state.vel[0]= rotvel
>
>
>
> #############################################################################
> # Run the simulation
>
> #############################################################################
>
>
> O.dt=0.5*PWaveTimeStep()
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([
> Bo1_Sphere_Aabb()
> ,Bo1_Box_Aabb()
> ], allowBiggerThanPeriod=True
> ),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom6D()
> ,Ig2_Box_Sphere_ScGeom6D()
> ],
> [Ip2_FrictMat_FrictMat_FrictPhys()
>
> ,Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,label="cohesiveIp")
> ],
> [Law2_ScGeom_FrictPhys_CundallStrack()
> ,Law2_ScGeom6D_CohFrictPhys_CohesionMoment()
> ]
> ), # End InteractionLoop
> NewtonIntegrator(damping=0.8,gravity=(0.,0.,0.)),
>
> ]
>
>
>
>
> --
> You received this question notification because your team yade-users 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
>

VG (varun-gupta) said : #2

Hello Jan,

Thanks for your response. The script seems to be running fine. Is there a problem with the interactions in my model set up ?

Regards
Varun

Jan Stránský (honzik) said : #3

Hi Varun,
the script works fine from the point of view that there are no errors.

However, in your question you describe cohesive bonds. If I got your
question correctly, your fragment is a set of particles mutually connected
with cohesive bonds. After running your script, there are no cohesive
bonds. This was my point.
Thanks
Jan

2016-05-02 1:22 GMT+02:00 VG <email address hidden>:

> Question #292911 on Yade changed:
> https://answers.launchpad.net/yade/+question/292911
>
> VG posted a new comment:
> Hello Jan,
>
> Thanks for your response. The script seems to be running fine. Is there
> a problem with the interactions in my model set up ?
>
> Regards
> Varun
>
> --
> You received this question notification because your team yade-users 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
>

VG (varun-gupta) said : #4

Hello Jan,

I have defined sample_material as CohFrictMat and also using Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,label="cohesiveIp") in the interaction loop. I am not sure what am I missing here. Could you please help me identify the problem in my model set up.

Thanks
Varun

Jan Stránský (honzik) said : #5

Hi Varun,
sorry, I was just stupid, after O.step() the interactions are there :-)

adding following functions at the end of the script (after O.step() )
should do what you want. You can use it also during simulation. Note that
this approach is probably not very efficient, but should work. In case the
inefficiency is too much, a more sophisticated graph theory should give
better answer.

In my case, I got several "aggregates" composed of 4 particles, but always
only 3 of them were connected by cohesively and one was standalone, which
corresponded to the result of aggregs() function.

cheers
Jan

PS: used approach is a nice exercise of recursive function call :-)

#########################################
def addBodyToAggreg(body,aggreg): # auxiliary function, add body [yade.Body
instance] and all its neighbors into aggreg (python set instance)
if body.id in aggreg: # do nothing if b is already in aggreg ...
return
aggreg.add(body.id) # ... otherwise add it to aggreg
intrs = body.intrs()
for i in intrs: # and add also all its neighbors ...
if not isinstance(i.phys,CohFrictPhys): # ... but only that connected with
CohFrictPhys interactions
continue
if i.phys.cohesionBroken: # ... but only not yet broken
continue
i2 = i.id1 if i.id2==body.id else i.id2 # choose the other body of
interaction
b2 = O.bodies[i2]
addBodyToAggreg(b2,aggreg) # and add it and all its neighbors to aggreg

def aggregs(): # actual function to detect standalone aggregates
ids = set(b.id for b in O.bodies if isinstance(b.shape,Sphere)) # first
make a set of all spheres ids
ret = []
while len(ids)>0: # while there are still some particles not assigned to
any aggregate ...
i = ids.pop() # ... choose one random ...
b = O.bodies[i]
a = set() # ... create new aggregate (set of sphere ids)
addBodyToAggreg(b,a) # ... and add the sphere together with all its
neigbors to aggregate
for bid in a: # delete all used ids from ids
ids.discard(bid)
ret.append(a)
return ret

aggs = aggregs()
for a in aggs:
print a
#########################################

2016-05-02 18:27 GMT+02:00 VG <email address hidden>:

> Question #292911 on Yade changed:
> https://answers.launchpad.net/yade/+question/292911
>
> VG posted a new comment:
> Hello Jan,
>
> I have defined sample_material as CohFrictMat and also using
>
> Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,label="cohesiveIp")
> in the interaction loop. I am not sure what am I missing here. Could you
> please help me identify the problem in my model set up.
>
> Thanks
> Varun
>
> --
> You received this question notification because your team yade-users 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
>

VG (varun-gupta) said : #6

Hello Jan,

Thanks a lot for providing these functions!

I am noticing something strange though, may be I am misinterpreting it. Running the same script with a seed value of 1 in makeCloud, leads to 19 initial aggregates with 11 particles each. I sent you the database file generated using memoizeDb in randomDensePack through email.

Now, the above functions print all the standalone aggregates. I am noticing that some of these aggregates are showing more than 11 particles. This shouldn't happen unless there is formation of new cohesive bonds. I am unable to explain this behavior.

Regards
Varun

Can you help with this problem?

Provide an answer of your own, or ask VG for more information if necessary.

To post a message you must log in.