apparent instability

This is a continuation of my previous question, "problem with bonds breaking".

I installed ESys version 2.3, and the bond instability seems to still be there. The bond stength is very high; bonds should not break in this simulation. The time step has been set 10 times lower than the (estimated) critical time step.

The code is:

#import the appropriate ESyS-Particle modules:
import sys, math
from numpy import *
from scipy import *
from esys.lsm import *
from esys.lsm.util import *
from esys.lsm.geometry import *
import random
from POVsnaps import POVsnaps
from POVsnaps2 import POVsnaps2
from POVsnaps3 import POVsnaps3
from WallLoader import WallLoaderRunnable

# instantiate a simulation object
# and initialise the neighbour search algorithm:
sim = LsmMpi(numWorkerProcesses=1, mpiDimList=[1,1,1])
sim.initNeighbourSearch(particleType="RotSphere",gridSpacing=2.5,verletDist=0.5)

# set values of simulation variables. Units are mks
numIceParticlesx = 40
numIceParticlesz = 12
particleRad = 0.5
particleRho = 900.0
particleMass = 4*pi*particleRad**3*particleRho/3.0 # for 3D
normalStiffness = 9000000.0 # reduced by a factor of 1000 to increase timestep
poisRatio = 0.25
kinFriction = 0.4
statFriction = 0.6
shearStrength = 5000000
tanPhi = 0.577

wallposition = -2.0*numIceParticlesx*particleRad - 0.6
print 'wall position =', wallposition

numsteps = 3200000
#stepsize = 0.2*particleRad*sqrt(particleRho/normalStiffness)
timestep = 0.0001
timeincrement = 32000

#set the number of timesteps and timestep increment:
sim.setNumTimeSteps(numsteps)
sim.setTimeStepSize(timestep)

#specify the spatial domain for the simulation:
domain = BoundingBox(Vec3(-120,-120,-120), Vec3(120,120,120))
sim.setSpatialDomain(domain)

# add the particles to the domain:
particleList = []
for countz in range(numIceParticlesz):
   zDist = -1.0*countz
   for countx in range(numIceParticlesx):
      xDist = -1.0*countx
      yDist = 0.0
# yDist = 0.001*(5-random.randint(1,10))
      particle_number = countz*100 + countx
      particle=SimpleSphere(id=particle_number, posn=Vec3(xDist,yDist,zDist), radius=particleRad, mass=particleMass)
# particle.setLinearVelocity(Vec3(0.0,0.0,0.0))
      sim.createParticle(particle)
      particleList.append(particle)

# create bonds between particles separated by less than the specified maxDist
sim.createConnections(ConnectionFinder(maxDist=0.005, bondTag=1, pList=particleList))

# create rotational elastic-brittle bonds between particles
pp_bonds = sim.createInteractionGroup(BrittleBeamPrms(name = "pp_bonds", youngsModulus=normalStiffness, poissonsRatio=poisRatio, cohesion=shearStrength, tanAngle=tanPhi, tag=1))

# initialize frictional interactions for unbonded particles
sim.createInteractionGroup(FrictionPrms(name="friction", youngsModulus=normalStiffness, poissonsRatio=poisRatio, dynamicMu=kinFriction, staticMu=statFriction))

# create an exclusion between bonded and frictional interactions
sim.createExclusion(interactionName1="pp_bonds", interactionName2="friction")

# add a vertical wall to push the particles along
sim.createWall(name = "icewall", posn=Vec3(wallposition,0,0), normal=Vec3(1,0,0))

# specify the type of interactions between wall and particles:
sim.createInteractionGroup(NRotElasticWallPrms(name = "elasticWall1", wallName = "icewall", normalK = normalStiffness))

# create mesh walls: ramp
sim.readMesh(fileName = "ramp.msh", meshName = "ramp_mesh_wall")

# specify that the particles undergo repulsion with the mesh wall
sim.createInteractionGroup(NRotElasticTriMeshPrms(name = "mesh_repel", meshName = "ramp_mesh_wall", normalK = normalStiffness))

# add a wall loader to move the icewall
wall_loader1 = WallLoaderRunnable (LsmMpi = sim, wallName = "icewall", vPlate = Vec3(0.1, 0.0, 0.0), startTime = 0, rampTime = 50)
sim.addPreTimeStepRunnable (wall_loader1)

#initialise gravity in the domain:
sim.createInteractionGroup(GravityPrms(name="earth-gravity", acceleration=Vec3(0,-9.81,0)))

# add buoyancy
sim.createInteractionGroup(BuoyancyPrms(name="water", acceleration=Vec3(0,-9.81,0), fluidDensity=1000, fluidHeight = 0.0))

#add local viscosity to simulate fluid resistance:
sim.createInteractionGroup(LinDampingPrms(name="linDamping",viscosity=0.95,maxIterations=100))

# add a checkpointer
sim.createCheckPointer(CheckPointPrms(fileNamePrefix = 'snapshot', beginTimeStep = 0, endTimeStep = numsteps, timeStepIncr = timeincrement))

#add a POVsnaps Runnable:
povcam = POVsnaps(sim=sim, interval=timeincrement)
#povcam.configure(lookAt=Vec3(0,0,-25), camPosn=Vec3(20,0,-25), zoomFactor=0.2)
povcam.configure(lookAt=Vec3(-20,0,-10), camPosn=Vec3(-20,10,-10), zoomFactor=0.02)
sim.addPostTimeStepRunnable(povcam)

#add another POVsnaps Runnable:
povcam2 = POVsnaps2(sim=sim, interval=timeincrement)
povcam2.configure(lookAt=Vec3(0,0,-2), camPosn=Vec3(10,0,-2), zoomFactor=0.08)
sim.addPostTimeStepRunnable(povcam2)

#add another POVsnaps Runnable:
povcam3 = POVsnaps3(sim=sim, interval=timeincrement)
povcam3.configure(lookAt=Vec3(0,0,0), camPosn=Vec3(0,0,15), zoomFactor=0.06)
sim.addPostTimeStepRunnable(povcam3)

#execute the simulation
sim.run()

The mesh file is:

Triangle
3D-Nodes 4
0 0 0 1.0 -0.5 10.0
1 1 0 1.0 -0.5 -10.0
2 2 0 61.0 29.5 10.0
3 3 0 61.0 29.5 -10.0

Tri3 2
0 0 1 3 0
1 0 0 3 2