Particle size distribution curve after particle breaking

Asked by yangjunjie on 2018-11-05

Hi ,

This is my last question about particle crushing https://answers.launchpad.net/yade/+question/675477.

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:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
2018-11-07
Last query:
2018-11-07
Last reply:
2018-11-06
Robert Caulk (rcaulk) said : #1

Hello,

I guess you want something like this:

def getRadiusArray():
    radii=[]
    for b in O.bodies:
 if isinstance(b.shape,Sphere):
         radii.append(b.shape.radius)

getRadiusArray()
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 : #2

**noticed it was missing some stuff**

def getRadiusArray():
....

    return radii

radii = getRadiusArray()
values,base=np.histogram(radii,bins=40)

yangjunjie (yangjunjie11) said : #3

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:

readParamsFromTable(rMean=.05,rRelFuzz=.3,maxLoad=1e6,minLoad=1e4)
from yade.params.table import *

from yade import pack, plot
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
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'),
   #PyRunner(command='getRadiusArray()',iterPeriod=100,label ='psd'),

]
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
    r = b.shape.radius
    # 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
#def getRadiusArray():
# radii=[]
# for b in O.bodies:
# if isinstance(b.shape,Sphere):
# radii.append(b.shape.radius)
# return radii
#
#radii=getRadiusArray()
#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.append(wall(max([b.state.pos[2]+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)
   O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]
   checker.command='unloadPlate()'

def unloadPlate():
   if abs(O.forces.f(plate.id)[2])>maxLoad:
      plate.state.vel*=-1
      checker.command='stopUnloading()'

def stopUnloading():
   if abs(O.forces.f(plate.id)[2])<minLoad:
      plot.saveDataTxt(O.tags['d.id']+'.txt')
      O.pause()

def addPlotData():
   if not isinstance(O.bodies[-1].shape,Wall):
      plot.addData(); return
   Fz=O.forces.f(plate.id)[2]
   plot.addData(Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],unbalanced=unbalancedForce(),i=O.iter)

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

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

waitIfBatch()

Best Robert Caulk (rcaulk) said : #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:
> https://answers.launchpad.net/yade/+question/675926
>
> 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:
>
> readParamsFromTable(rMean=.05,rRelFuzz=.3,maxLoad=1e6,minLoad=1e4)
> from yade.params.table import *
>
> from yade import pack, plot
> O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
> 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'),
> #PyRunner(command='getRadiusArray()',iterPeriod=100,label ='psd'),
>
> ]
> 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
> r = b.shape.radius
> # 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
> #def getRadiusArray():
> # radii=[]
> # for b in O.bodies:
> # if isinstance(b.shape,Sphere):
> # radii.append(b.shape.radius)
> # return radii
> #
> #radii=getRadiusArray()
> #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.append(wall(max([b.state.pos[2]+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)
> O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]
> checker.command='unloadPlate()'
>
> def unloadPlate():
> if abs(O.forces.f(plate.id)[2])>maxLoad:
> plate.state.vel*=-1
> checker.command='stopUnloading()'
>
> def stopUnloading():
> if abs(O.forces.f(plate.id)[2])<minLoad:
> plot.saveDataTxt(O.tags['d.id']+'.txt')
> O.pause()
>
> def addPlotData():
> if not isinstance(O.bodies[-1].shape,Wall):
> plot.addData(); return
> Fz=O.forces.f(plate.id)[2]
>
> plot.addData(Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],unbalanced=unbalancedForce(),i=O.iter)
>
> plot.plots={'i':('unbalanced',),'w':('Fz',)}
> plot.plot()
>
> O.run()
> #from yade import utils
> #binsSizes, binsProc, binsSumCum= utils.psd(bins=5, mass=True)
>
> waitIfBatch()
>
> --
> You received this question notification because you are subscribed to
> the question.
>

yangjunjie (yangjunjie11) said : #5

Thanks Robert Caulk, that solved my question.

Robert Caulk (rcaulk) said : #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:

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:

readParamsFromTable(rMean=.05,rRelFuzz=.3,maxLoad=1e6,minLoad=1e4)
from yade.params.table import *

from yade import pack, plot
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
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'),
   #PyRunner(command='getRadiusArray()',iterPeriod=100,label ='psd'),

]
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
    r = b.shape.radius
    # 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
#def getRadiusArray():
# radii=[]
# for b in O.bodies:
# if isinstance(b.shape,Sphere):
# radii.append(b.shape.radius)
# return radii
#
#radii=getRadiusArray()
#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.append(wall(max([b.state.pos[2]+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)
   O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=200)]
   checker.command='unloadPlate()'

def unloadPlate():
   if abs(O.forces.f(plate.id)[2])>maxLoad:
      plate.state.vel*=-1
      checker.command='stopUnloading()'

def stopUnloading():
   if abs(O.forces.f(plate.id)[2])<minLoad:
      plot.saveDataTxt(O.tags['d.id']+'.txt')
      O.pause()

def addPlotData():
   if not isinstance(O.bodies[-1].shape,Wall):
      plot.addData(); return
   Fz=O.forces.f(plate.id)[2]
   plot.addData(Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],unbalanced=unbalancedForce(),i=O.iter)

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

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

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

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 : #8

Sorry, that's my fault.

Robert Caulk (rcaulk) said : #9

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

[1]https://yade-dem.org/doc/yade.utils.html?highlight=utils.psd#yade.utils.psd