Compression of a packing created using RegularOrtho

Hi,

I performed a uniaxial compression of a packing created using regularOrtho (https://yade-dem.org/doc/yade.pack.html?highlight=pack#yade.pack.regularOrtho). But this boned particle model explodes when it was compressed at a certain time. Honestly, I intend to simulate crack initiation and propagation of BPM under compression loads.
Can you tell me what happed to my model ? Here is a sgement of my code:
Thanks in advance.
Xavier

####### My code #######
# !/usr/bin/python
# -*- coding: utf-8 -*-

from yade import utils
from yade import pack
from yade import plot
from yade import qt

# input parameters
porosity = 0.4764012244
Young = 1.0e11 # 100Gpa
FrictAng = radians(18)
Density = 2640/(1-porosity)
Poisson = 1/3.0
Cohesion = 4.5e7 # pa
TensileStr = 4.5e6 # pa

rock = JCFpmMat(young=Young,frictionAngle=FrictAng,density=Density,poisson=Poisson,tensileStrength=TensileStr,cohesion=Cohesion)
steel = JCFpmMat(young=2.06e11,frictionAngle=radians(30),density=7800,poisson=0.3,tensileStrength=0,cohesion=0)
for Mat in (rock,steel):
 O.materials.append(Mat)

# regularOrth packing
pred = pack.inAlignedBox((0,0,0),(0.01,0.03,0.06))
assembly = pack.regularOrtho(pred,radius=0.001,gap=0,material=rock)
O.bodies.append(assembly)

# aabbWalls
wall = O.bodies.append(utils.aabbWalls(material=steel,thickness=0.001,color=(1,1,1)))

# delete other walls only save walls in z axis
O.bodies.erase(wall[0])
O.bodies.erase(wall[1])
O.bodies.erase(wall[2])
O.bodies.erase(wall[3])

# fixed the bottom wall
bottom_wall = O.bodies[wall[4]]
bottom_wall.state.blockedDOFs = 'xyzXYZ'

# set velocity for upper wall
upper_wall = O.bodies[wall[5]]
upper_wall.state.blockedDOFs = 'xyzXYZ'
upper_wall.state.vel = (0,0,-0.005) # 5mm/s

###########################
# ENGINES
###########################

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.25),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.25),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,label='lawFunctor')]
 ),
 # VTKRecorder(fileName='post/uniCom-',recorders=['all'],iterPeriod=50),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=5,timestepSafetyCoefficient=0.8,defaultDt=PWaveTimeStep()),
 NewtonIntegrator(damping=0.3,gravity=(0,0,-9.81)),
]
O.trackEnergy=True
O.step()

#set an optimal timestep
O.dt = utils.PWaveTimeStep()
O.usesTimeStepper = True

#3D view and controller
try:
 qt.Controller()
 qt.View()
except:
 pass