create a sample with no voids between particles,so it seems like a solid

Hello to everyone,

i have this code, and i want to simulate the cutting process of a rock sample.
Above is my code, which generate a sample with specific characteristics and a
cutter (geom.facetBox).

When i put my data in paraview, i see that are many gaps between spheres. It seems like my sample,
is more like a gas rather as a rock sample.

I want to simulate the fracture pattern and the chip formation during cutting, byt i cant!!!
Question: how can i create a sample with no voids? I try periodic triaxial test byt there are voids too.

My code is here:

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division

from yade import plot,pack,timing, utils, geom
import time, sys, os, copy

#import matplotlib
#matplotlib.rc('text',usetex=True)
#matplotlib.rc('text.latex',preamble=r'\usepackage{concrete}\usepackage{euler}')

# default parameters or from table
readParamsFromTable(noTableOk=True, # unknownOk=True,
 young=24e9, #kn
 poisson=.2, #ks

 sigmaT=3.5e6,
 frictionAngle=atan(0.8), #frictionAngle
 epsCrackOnset=1e-4,
 relDuctility=30,

 intRadius=1.5,
 dtSafety=.8,
 damping=0.4,
 specimenLength=1,
 sphereRadius=0.99e-3, #0.35e-3

 # isotropic confinement (should be negative)
 isoPrestress=0,
)

from yade.params.table import *

if 'description' in O.tags.keys(): O.tags['id']=O.tags['id']+O.tags['description']
#material properties
concreteId=O.materials.append(CpmMat(young=young,frictionAngle=frictionAngle,poisson=poisson,density=4800,sigmaT=sigmaT,
     relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress))

sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((-75e-3*specimenLength,-27e-3*specimenLength,-6.35e-3*specimenLength),(75e-3*specimenLength,27e-3*specimenLength,6.35e-3*specimenLength)),spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
sp.toSimulation(material=concreteId)

bx = geom.facetBox(center=(-80e-3,27e-3,0.0),extents=(3e-3,-5e-3,6.35e-3),orientation=Quaternion((0,0,1),-pi/36), wallMask=(2+8))
O.bodies.append(bx)
for facet in bx:
 facet.state.blockedDOFs='xyzXYZ'
 facet.state.vel=(5.0,0,0)
#bx.state.blockedDOFs = 'xyzXYZ'
#bx.state.vel = (10.0,0,0)

O.dt=dtSafety*PWaveTimeStep()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),Bo1_Facet_Aabb()],verletDist=.05*sphereRadius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc'), Ig2_Facet_Sphere_ScGeom()],
  [Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=1)],
  [Law2_ScGeom_CpmPhys_Cpm()],
 ),
 NewtonIntegrator(damping=damping,label='damper'),
    VTKRecorder(fileName='3d-vtk-',recorders=['cracks'],iterPeriod=1000),
    CpmStateUpdater(iterPeriod=50),
    PyRunner(command='addPlotData()',realPeriod=2),
    PyRunner(iterPeriod=100,initRun=True,command='recorder()'),
]

def addPlotData():
   global Fx
   Fx = 0.0
   Fx = abs(sum(O.forces.f(facet.id)[0] for facet in bx))
   plot.addData(i = O.iter , Fx=Fx)
   plot.saveDataTxt(O.tags['d.id']+'.txt')

tensCks=shearCks=cks=cks0=0
def recorder():
    global tensCks, shearCks
    tensCks=0
    shearCks=0
    for o in O.bodies:
        if isinstance(o.shape,Sphere):
            tensCks+=o.state.tensBreak
            shearCks+=o.state.shearBreak
    yade.plot.addData({ 't':O.time,'i':O.iter,'tc':tensCks,'sc':shearCks})
    plot.saveDataTxt('intrs.txt')

plot.plots={'i':('Fx')}
plot.plot()

Thank you!!