← 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 proposed the following answer:
This is a good python teaching moment. Good luck ;-)

Le mar. 6 nov. 2018 à 04:43, yangjunjie <
question675926@xxxxxxxxxxxxxxxxxxxxx> 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.
>

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