← Back to team overview

yade-users team mailing list archive

Re: [Question #675926]: Particle size distribution curve after particle breaking

 

Question #675926 on Yade changed:
https://answers.launchpad.net/yade/+question/675926

Robert Caulk posted a new comment:
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()
>>>>>>>>>>>>>>>>>>>>>

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.