Clumping - angl of repose

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 20000 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=10000)
# 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
   plot.addData(i=O.iter,unbalanced=unbalancedForce(),**O.energy)
#O.save('Modifiedd.txt.bz2')
#while 1:
 #O.run(100,True)
 #if unbalancedForce()<1e-5:
  #break
plot.plots={'i':('unbalanced',None,O.energy.keys)}
plot.plot()

O.saveTmp()

#plot.saveDataTxt('2')
from yade import qt
qt.View()