yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #18346
Re: [Question #675926]: Particle size distribution curve after particle breaking
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 your team yade-users is
an answer contact for Yade.