# Particle size distribution curve after particle breaking

Hi ,

I use fragment replace method to simulate particle breaking (the new generation of smaller particles that replace the disappearing broken particle), but I must to observe the evolving of particle size distribution and changing between initial PSD curve and ultimate one.
The x-axis of particle size distribution curve is particle size(Logarithmic coordinates) and y-axis is mass percentage less than a certain particle size.
So how can I plot this curve ?

Thanks

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
2018-11-07
Last query:
2018-11-07
2018-11-06
 Robert Caulk (rcaulk) said on 2018-11-05: #1

Hello,

I guess you want something like this:

for b in O.bodies:
if isinstance(b.shape,Sphere):

values, base = np.histogram(data, bins=40)
cumulative = np.array(np.cumsum(values),dtype=float)
perc = cumulative/max(cumulative)*100
plt.semilogx(base[:-1], cumulative, c='blue')

Cheers,

Robert

 Robert Caulk (rcaulk) said on 2018-11-05: #2

**noticed it was missing some stuff**

....

 yangjunjie (yangjunjie11) said on 2018-11-06: #3

Hi Robert:

I added the code you provided to mine, there seems to be some conflict between plt.show() and plot.plot(), Or my understanding is wrong. Help me plz..

This is my whole code:

sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=rMean)
sp.toSimulation()

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
PyRunner(command='replace()',iterPeriod=100,label ='rep'),

]
O.dt=.5*PWaveTimeStep()
##############################################################
def selectBodiesToCrush():
stresses = bodyStressTensors() # e.g. based on overall particle stress
ret = []
for i,stress in enumerate(stresses):
if stress[0,0] > 50: # Particle breakage criterion
ret.append(i)
return ret

def createChildren(b):
p = b.state.pos
# Particle replacement mode
s1 = sphere(p+(.5*r,0,0),.5*r)
s2 = sphere(p-(.5*r,0,0),.5*r)
return s1,s2

def replace():
bodiesToDelete = selectBodiesToCrush()
for i in bodiesToDelete:
if isinstance(O.bodies[i].shape,Sphere): # Limitation of the smallest particles that can break
b = O.bodies[i]
children = createChildren(b)
O.bodies.erase(i)
O.bodies.append(children)

##############################################################
#import matplotlib.pyplot as plt
#import numpy as np
# for b in O.bodies:
# if isinstance(b.shape,Sphere):
#
#cumulative = np.array(np.cumsum(values),dtype=float)
#perc = cumulative/max(cumulative)*100
#plt.semilogx(base[:-1], cumulative, c='blue')
#plt.show()
##############################################################

def checkUnbalanced():
if O.iter<5000: return
if unbalancedForce()>.1: return
O.bodies.append(wall(max([b.state.pos+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))
global plate
plate=O.bodies[-1]
plate.state.vel=(0,0,-.1)

plate.state.vel*=-1

plot.saveDataTxt(O.tags['d.id']+'.txt')
O.pause()

if not isinstance(O.bodies[-1].shape,Wall):
Fz=O.forces.f(plate.id)

plot.plots={'i':('unbalanced',),'w':('Fz',)}
plot.plot()

O.run()
#binsSizes, binsProc, binsSumCum= utils.psd(bins=5, mass=True)

waitIfBatch() Robert Caulk (rcaulk) said on 2018-11-06: #4

This is a good python teaching moment. Good luck ;-)

Le mar. 6 nov. 2018 à 04:43, yangjunjie <
<email address hidden>> a écrit :

> Question #675926 on Yade changed:
>
> yangjunjie posted a new comment:
> Hi Robert:
>
> Thank you for your reply !! But I meet some problems:
> I added the code you provided to mine, there seems to be some conflict
> between plt.show() and plot.plot(), Or my understanding is wrong. Help me
> plz..
>
> This is my whole code:
>
>
> from yade import pack, plot
> sp=pack.SpherePack()
> sp.makeCloud((0,0,0),(1,1,1),rMean=rMean)
> sp.toSimulation()
>
> O.engines=[
> ForceResetter(),
>
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
> InteractionLoop(
>
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
> [Ip2_FrictMat_FrictMat_FrictPhys()],
> [Law2_ScGeom_FrictPhys_CundallStrack()]
> ),
> NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
> PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
> PyRunner(command='replace()',iterPeriod=100,label ='rep'),
>
> ]
> O.dt=.5*PWaveTimeStep()
> ##############################################################
> def selectBodiesToCrush():
> stresses = bodyStressTensors() # e.g. based on overall particle stress
> ret = []
> for i,stress in enumerate(stresses):
> if stress[0,0] > 50: # Particle breakage criterion
> ret.append(i)
> return ret
>
> def createChildren(b):
> p = b.state.pos
> # Particle replacement mode
> s1 = sphere(p+(.5*r,0,0),.5*r)
> s2 = sphere(p-(.5*r,0,0),.5*r)
> return s1,s2
>
> def replace():
> bodiesToDelete = selectBodiesToCrush()
> for i in bodiesToDelete:
> if isinstance(O.bodies[i].shape,Sphere): # Limitation of the
> smallest particles that can break
> b = O.bodies[i]
> children = createChildren(b)
> O.bodies.erase(i)
> O.bodies.append(children)
>
> ##############################################################
> #import matplotlib.pyplot as plt
> #import numpy as np
> # for b in O.bodies:
> # if isinstance(b.shape,Sphere):
> #
> #values, base = np.histogram(radii, bins=40)
> #cumulative = np.array(np.cumsum(values),dtype=float)
> #perc = cumulative/max(cumulative)*100
> #plt.semilogx(base[:-1], cumulative, c='blue')
> #plt.show()
> ##############################################################
>
> def checkUnbalanced():
> if O.iter<5000: return
> if unbalancedForce()>.1: return
> O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))
> global plate
> plate=O.bodies[-1]
> plate.state.vel=(0,0,-.1)
>
> plate.state.vel*=-1
>
> plot.saveDataTxt(O.tags['d.id']+'.txt')
> O.pause()
>
> if not isinstance(O.bodies[-1].shape,Wall):
> Fz=O.forces.f(plate.id)
>
>
> plot.plots={'i':('unbalanced',),'w':('Fz',)}
> plot.plot()
>
> O.run()
> #binsSizes, binsProc, binsSumCum= utils.psd(bins=5, mass=True)
>
> waitIfBatch()
>
> --
> the question.
>

 yangjunjie (yangjunjie11) said on 2018-11-07: #5

Thanks Robert Caulk, that solved my question.

 Robert Caulk (rcaulk) said on 2018-11-07: #6

I'm sorry, I disagree with your choice to edit the thread after you received help from the community. It may come as a surprise to you, but others will want to use this thread in the future. Therefore, I am going to go ahead and re-post the response that you removed:

>>>>>>>>>>>>>>>>>>>>>
Hi Robert:

I added the code you provided to mine, there seems to be some conflict between plt.show() and plot.plot(), Or my understanding is wrong. Help me plz..

This is my whole code:

sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=rMean)
sp.toSimulation()

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
PyRunner(command='replace()',iterPeriod=100,label ='rep'),

]
O.dt=.5*PWaveTimeStep()
##############################################################
def selectBodiesToCrush():
stresses = bodyStressTensors() # e.g. based on overall particle stress
ret = []
for i,stress in enumerate(stresses):
if stress[0,0] > 50: # Particle breakage criterion
ret.append(i)
return ret

def createChildren(b):
p = b.state.pos
# Particle replacement mode
s1 = sphere(p+(.5*r,0,0),.5*r)
s2 = sphere(p-(.5*r,0,0),.5*r)
return s1,s2

def replace():
bodiesToDelete = selectBodiesToCrush()
for i in bodiesToDelete:
if isinstance(O.bodies[i].shape,Sphere): # Limitation of the smallest particles that can break
b = O.bodies[i]
children = createChildren(b)
O.bodies.erase(i)
O.bodies.append(children)

##############################################################
#import matplotlib.pyplot as plt
#import numpy as np
# for b in O.bodies:
# if isinstance(b.shape,Sphere):
#
#cumulative = np.array(np.cumsum(values),dtype=float)
#perc = cumulative/max(cumulative)*100
#plt.semilogx(base[:-1], cumulative, c='blue')
#plt.show()
##############################################################

def checkUnbalanced():
if O.iter<5000: return
if unbalancedForce()>.1: return
O.bodies.append(wall(max([b.state.pos+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1))
global plate
plate=O.bodies[-1]
plate.state.vel=(0,0,-.1)

plate.state.vel*=-1

plot.saveDataTxt(O.tags['d.id']+'.txt')
O.pause()

if not isinstance(O.bodies[-1].shape,Wall):
Fz=O.forces.f(plate.id)

plot.plots={'i':('unbalanced',),'w':('Fz',)}
plot.plot()

O.run()
#binsSizes, binsProc, binsSumCum= utils.psd(bins=5, mass=True)

waitIfBatch()
>>>>>>>>>>>>>>>>>>>>>

 Bruno Chareyre (bruno-chareyre) said on 2018-11-07: #7

Yes, please, yangjunjie, don't re-write the history of discussions.
I'm also curious why the first part of the question was removed: "This is my last question about particle crushing".
It seems last question on crushing ends up as first question on the additivity of python codes. ;)

Please think twice _before_ posting, not after posting.

Bruno

 yangjunjie (yangjunjie11) said on 2018-11-08: #8

Sorry, that's my fault.

 Robert Caulk (rcaulk) said on 2018-11-08: #9

Turns out, Yade might have a nice built in function for this. I haven't used it myself, but it could be worth trying.