Errors during step-by-step particle deletion of consolidated specimen

Hi all,

I am simulating the step-by-step fine particle deletion process of a consolidated gap-graded specimen in 2D case. The confining pressure is maintained constant during the deletion process. The initial fine content by mass is FC = 25%, and the 1% of fine particles are randomly deleted at each loop. After 25 loops, all fine particles should be deleted as expected.

At the first several loop of my case, fine particles can be deleted correctly. However, the simulation suspend after fine loss = 8% though there still many fines in the specimen. It seems that the confining pressure cannot be re-achieved after several particle deletions? Attached please see a simple case to reproduce the ERROR. Could you please help to see where is the problem.

Many thanks and best regards,
Zheng

#### Code below ####

from yade import pack,plot,export
import matplotlib.pyplot as plt
import numpy as np
import random

O.materials.append(FrictMat(young=6.e8,poisson=.8,frictionAngle=.0))

sigmaIso=-1e5
sp = pack.SpherePack()
size = 0.24
# 25% fines in 2d case.
sp.makeCloud(minCorner=(0,0,.05),maxCorner=(size,size,.05),psdSizes=[0.001,0.0012,0.005,0.006], psdCumm=[0,0.0625,0.0625,1.0],periodic=True,seed=1,distributeMass=True,num=1500)

maxFine = 0.0012 # maximum diameter of fine particles
sp.toSimulation()
O.cell.hSize = Matrix3(size,0,0, 0,size,0, 0,0,.1) # used for periodic boundaries.
massAll = 0
massEroded = 0

for p in O.bodies:
   p.state.blockedDOFs = 'zXY'
   p.state.mass = 2650 * 0.1 * pi * p.shape.radius**2 # 0.1 = thickness
   inertia = 0.5 * p.state.mass * p.shape.radius**2
   p.state.inertia = (.5*inertia,.5*inertia,inertia)
   massAll += p.state.mass # mass of all particles

O.dt = utils.PWaveTimeStep()

O.engines = [
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   PeriTriaxController(
      dynCell=True,
      goal=(-1e5,-1e5,0),
      stressMask=3,
      relStressTol=.001,
      maxUnbalanced=.1,
      maxStrainRate=(.5,.5,.0),
      doneHook='delFine()',
      label='biax'
   ),
   NewtonIntegrator(damping=.1),
]

delperc = range(1,25,1)
delperc = [x/100 for x in delperc]
n = 0
def delFine():
   global delperc,n,massEroded
   #global n
   #global massEroded
   if n < len(delperc):
      print('Current stress before deletion, ',getStress())
      setContactFriction(0.5)
      bodyRadius=[]
      for b in O.bodies:
        if b.shape.radius<=maxFine/2:
         bodyRadius.append([b.id,b.shape.radius,b.state.mass])
      bodyRadius.sort(key=lambda x:x[1])
      i = 0
      for b in bodyRadius:
        if massEroded <= delperc[n]*massAll:
           global massEroded
           massEroded += b[2]
           O.bodies.erase(b[0])
           i+=1
        else:
           break
      print('The amount of fines loss by number, ',i)
      print('Current stress after deletion = ',getStress())
      # when the simulation finish?
      if i == 0:
        print(delperc[n]*100,'% fines loss!')
        O.save(str(int(delperc[n]*100))+'.yade.gz')
        n += 1
   if n == len(delperc):
      biax.doneHook='Finished()'

def Finished():
   print('Test Finished')
   O.pause()

O.run()
O.wait()