New question #291202 on Yade:
https://answers.launchpad.net/yade/+question/291202
Hi all,
In below code I have tried to study the effect of percentage/ratio of clumping
on repose angle. before adding the clumping section the scrip worked properly
however by adding clump grain do not fall down because of their gravity and the
sample explodes .
Would you please advise?
Thanks
#!/usr/bin/python
from yade import pack,utils, qt, plot
pred = pack.inAlignedBox((0,0,0),(14.5,14.5,14.5))
#create material
soil1 =
CohFrictMat(young=30e9,poisson=0.3,frictionAngle=radians(30),density=2600.0,normalCohesion=1e6,
shearCohesion=80e6,label='soil')
#color=(1,0,0) red color
#soil1 =
FrictMat(young=1e6,poisson=0.4,frictionAngle=radians(30),density=2500.0,label='soil')
O.materials.append(soil1)
O.bodies.append(utils.wall(0,axis=1,sense=1))
O.materials.append(CohFrictMat(young=30e9,poisson=0.3, frictionAngle =
radians(30) , label='wallmat'))
wallmat = O.materials[-1]
45
spheres=SpherePack()
spheres=pack.randomDensePack(pred,radius=.5,material='soil',spheresInCell=1000,color=(1,0,0),returnSpherePack=True)
spheres.toSimulation()
#O.bodies.append(spheres)
#
#
O.engines=[
ForceResetter(),#reset forces
InsertionSortCollider([Bo1_Wall_Aabb(),Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Wall_Sphere_ScGeom()], # collision geometry
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()], #
collision "physics"
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()]
# contact law -- apply forces
),
# apply gravity force to particles
# damping: numerical dissipation of energy
NewtonIntegrator(damping=0.5,gravity=(0,-9.81,0)),
#qt.SnapshotEngine(fileBase='3d-',iterPeriod=200,label='snapshot'),
# this engine will be called after 2 steps, only once
PyRunner(command='checkUnbalanced()',realPeriod=200,label='checker'),
PyRunner(command='addPlotData()',iterPeriod=200)
]
#CLUMPING###
#create a box:
#id_box =
O.bodies.append(box((0,0,0),(.0045,0.0045,.0045),fixed=True,material=Mat))
#create assembly of spheres:
#sp=pack.SpherePack()
#sp.makeCloud((0,0,0),(.004,0.004,.004),rMean=.00018,rRelFuzz=.0001,num=1000,periodic=False)
#O.bodies.append([sphere(c,r,material=Mat) for c,r in sp])
#print len(sp),' particles generated.'
print 'Roundness coefficient without clumps is: ',O.bodies.getRoundness()
show how to use makeClumpTemplate():
#dyad:
#relRadList1 = [.0002,.0002]
#relPosList1 = [[.0003,0,0],[0.0006,0,0]]
#peanut:
relRadList2 = [.0002,0.0004,0.0002]
relPosList2 = [[0.0003,0,0],[0,0,0],[0.0006,0,0]]
#stick:
#relRadList3 = [1,1,1,1,1]
#relPosList3 = [[0,1,0],[0,2,0],[0,3,0],[0,4,0],[0,5,0]]
templates= []
#templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1))
templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
#templates.append(clumpTemplate(relRadii=relRadList3,relPositions=relPosList3))
show how to use replaceByClumps():
#replace by 50% dyads, 30% peanuts and 10% sticks:
#O.bodies.replaceByClumps(templates,[.5,.3,.1])
O.bodies.replaceByClumps(templates,[0.5])
show how to use getRoundness():
#create a list of all standalone spheres:
standaloneList = []
for b in O.bodies:
if b.isStandalone:
standaloneList.append(b.id)
print 'Roundness coefficient for spheres and clumps is:
',O.bodies.getRoundness()
print 'Roundness coefficient just for clumps is:
',O.bodies.getRoundness(standaloneList)
#sp.toSimulation()
print 'Number of elements: ', len(O.bodies)
print 'Box Volume: '
print 'Clumprepose,0.5,.25,30,3grain,0.5%'
O.trackEnergy=True
# set timestep to a fraction of the critical timestep
# the fraction is very small, so that the simulation is not too fast
# and the motion can be observed
O.dt=1*utils.PWaveTimeStep()
#makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=1)
# save the simulation, so that it can be reloaded later, for experimentation
#O.saveTmp()0
#checker.command='stopUnloading()'
def checkUnbalanced():
if unbalancedForce()<5e-2:
print('Reached target , stopping')
O.pause()
plot.saveDataTxt('Clumprepose,0.5,.5,30,3grain,0.5%.bz2')
plot.saveDataTxt('Clumprepose,0.5,.5,30,3grain,0.5%.txt')
#plot.saveDataTxt('modi.data.bz2')
# plot.saveGnuplot('bbb') is also possible
# collect history of data which will be plotted
def addPlotData():
# each item is given a names, by which it can be the unsed in plot.plots
# the **O.energy converts dictionary-like O.energy to plot.addData arguments
