How to glue spheres to the box?

Dear all,

I am trying to simulate direct shear experiments.

In my simulation, the periodic boundary condition is applied and in addition, the top and bottom are confined by a box.
To simulate a rough boundary surface similar to the experiment conditions, I tried to glue some spheres on the box.
I thought "clump" spheres and boxes together could be a feasible way. However, it seems not to work (or maybe I did it wrong).
It turns out the spheres penetrate into the box during isotropic compression.
I also tried to glue them together by playing with cohesive, but the contact model I use is a non-cohesive Hertz contact model.

I am wondering if there are other ways I can glue the spheres and box together and therefore, the glued material can also act as confinement. Thanks!

Hereby I attach the simplified code:

##########################
from yade import pack
import math

sp=pack.SpherePack()

O.periodic=True

RADIUS=0.025
length=0.75
height=0.375
width=0.375
thickness=RADIUS
PI=1.e5
spFRIC=26.6

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

O.materials.append(FrictMat(density=2500,young=5.5e10,poisson=0.25,frictionAngle=radians(spFRIC),label='boxMat'))

upBoundary=[]
lowBoundary=[]

### Creat top and bottom box
upBoundary.append(O.bodies.append(utils.box(center=(length/2,2*height+thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')))
lowBoundary.append(O.bodies.append(utils.box(center=(length/2,height-thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,width),fixed=1,wire=True,color=(1,0,0),material='boxMat')))

### Roughness setting
wallLength=15
wallWidth=10
roughnessPara=2

### top spheres
for c in range(1, wallLength):
    for d in range(0, wallWidth):
        upBoundary.append(O.bodies.append(sphere([c*RADIUS*roughnessPara,2*height-thickness,d*RADIUS],RADIUS,fixed=1,color=(72/255.,72/255.,147/255.),material='boxMat')))

### low spheres
for a in range(1, wallLength):
    for b in range(0, wallWidth):
        lowBoundary.append(O.bodies.append(sphere([a*RADIUS*roughnessPara,height+thickness,b*RADIUS],RADIUS,fixed=1,color=(72/255.,72/255.,147/255.),material='boxMat')))

### Clump spheres and box together
idClump1=O.bodies.clump(upBoundary)
idClump2=O.bodies.clump(lowBoundary)

sp.makeCloud((0,height+3*RADIUS,0),(length,2*height-3*RADIUS,width), rMean=0.025, periodic=True, seed =1)
sphere_id = O.bodies.append([utils.sphere(s[0],s[1],color=(0,0,1),material='boxMat') for s in sp])

effCellVol=(O.bodies[idClump1].state.pos[1]-O.bodies[idClump2].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_FrictMat_FrictMat_MindlinPhys()],
  [Law2_ScGeom_MindlinPhys_Mindlin()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep())
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ]

def triaxDone():
 global phase
 vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
 contactStress=getStress(vol)
 vol_s=Rmean=Rmax=nbSph=0
 Rmin=1e6
 for o in O.bodies:
  if isinstance(o.shape,Sphere):
   nbSph+=1
   Rmean+=o.shape.radius
   if o.shape.radius>Rmax: Rmax=o.shape.radius
   if o.shape.radius<Rmin: Rmin=o.shape.radius
   vol_s += 4.*pi/3.*(o.shape.radius)**3
 Rmean=Rmean/nbSph
 O.pause()

#####################

Cheers,
Chien-Cheng