Decrease bond strengths to make a model collapse

Hello all,

I am trying to setup a simple simulation to evaluate the geometry of a slope after collapse. I have applied gravity to the model, let the particles settle and removed some elements to create the slope. What I want to do now is to decrease the bond strengths in order to make the slope collapse.

There is no any external force that makes the slope collapse, I just want to do it in an "artificial" way to evaluate the geometry. I am using a CohFrictMat and I am not sure how to do it. I have read in a paper that they do it by imposing lower tensile strengths to all bonds in the system, allowing bonds to break where the gravitationally induced stress exceeded the bond strength. Do I have to loop over O.interactions?

Thanks in advance.

Here is my code so far:

#!/usr/bin/python
from yade import pack,utils,plot
pred = pack.inAlignedBox((0,0,0),(14,10,8))
#create material
soil1 = CohFrictMat(young=1e6,poisson=0.4,frictionAngle=radians(30),density=2500.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=0))
O.bodies.append(utils.wall(8,axis=2,sense=0))
O.bodies.append(utils.wall(0,axis=2,sense=0))
O.bodies.append(utils.wall(14,axis=0,sense=0))
O.bodies.append(utils.wall(0,axis=0,sense=0))

spheres=SpherePack()
spheres=pack.randomDensePack(pred,radius=.5,material='soil',spheresInCell=200,color=(1,0,0),returnSpherePack=True)
spheres.toSimulation()

# simulation loop
O.engines=[
 ForceResetter(),#reset forces
           #collisions
           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.1,gravity=(0,-9.81,0)),
                    #check static equilibrium
                          PyRunner(command='checkUnbalanced()',realPeriod=2),
                    #plot unbalanced force vrs. iterations
                          PyRunner(command='addPlotData()',iterPeriod=100)
                        ]

def checkUnbalanced():
       if utils.unbalancedForce()<.05:
          O.pause()
          nx,ny,nz,d = 3.,2.,0.,-30 # for plane equation nx*x+ny*y+nz*z+d=0
          for b in O.bodies:
              p = b.state.pos
              if nx*p[0]+ny*p[1]+nz*p[2]+d>0 and p[1]>3: # p is on positive halfspace of plane nx*x+ny*y+n z*z+d=0
                 O.bodies.erase(b.id)

def addPlotData():
       # each item is given a names, by which it can be the unsed in plot.plot
        plot.addData(i=O.iter,unbalanced=utils.unbalancedForce())

# define how to plot data: 'i' (step number) on the x-axis, unbalanced force
# on the left y-axis, all energies on the right y-axis
# (O.energy.keys is function which will be called to get all defined energies)
# None separates left and right y-axis
plot.plots={'i':('unbalanced')}
# show the plot on the screen, and update while the simulation runs
plot.plot()

# 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=.5e-2*utils.PWaveTimeStep()

# save the simulation, so that it can be reloaded later, for experimentation
O.saveTmp()

from yade import qt
qt.View()
#O.run()