problem on biaxial compression test using clumps

Dear All,
I have a problem in simulating biaxial compression test using clumps. I have blocked the degrees of freedom as 'zXY' after initial isotropic consolidation. However, during shearing, it can be found that some particles still rotate along Y axis from the graphical view and the volumetric strain curve is also not correct. I wonder the reason.

ubuntu : 18.04
yade 2018.02b

There are minimum of my codes:

Code 1: preparation of clumps

from yade import pack
import random

sp = pack.SpherePack()
size = .32
sp.makeCloud(minCorner=(0,0,.05),maxCorner=(size,size,.05),rMean=.005,rRelFuzz=.2,num=400,periodic=True,seed=1)
sp.toSimulation()

ori = 0. * pi
coef = pi + 0.315 * sqrt(1-0.315**2) - acos(0.315)
fout = file('clump_zero.lst','w')

for b in O.bodies:
   pos = b.state.pos
   radius = b.shape.radius
   line_oo = 0.63 * radius
   eq_rad = sqrt(2 * coef / pi) * radius
   #print eq_rad
   x1 = pos[0] - 0.5 * line_oo * cos(ori)
   x2 = pos[0] + 0.5 * line_oo * cos(ori)
   y1 = pos[1] - 0.5 * line_oo * sin(ori)
   y2 = pos[1] + 0.5 * line_oo * sin(ori)
   fout.write(str(x1)+' '+str(y1)+' '+str(x2)+' '+str(y2)+' '+str(radius)+' '+str(eq_rad)+'\n')

fout.close()

Code 2: consolidation and shearing

from yade import export,plot,pack

mat_id = O.materials.append(FrictMat(young=4.e8,poisson=.8,frictionAngle=0.05,density=2670))
mat = O.materials[mat_id]

fin = open('clump_zero.lst','r')
clumpDataBase = fin.readlines()
fin.close()

z = 0.05
clump_id = []
coef = pi + 0.315 * sqrt(1-0.315**2) - acos(0.315)
area = 0.0
for i in xrange(len(clumpDataBase)):
   x1 = float(clumpDataBase[i].split()[0])
   y1 = float(clumpDataBase[i].split()[1])
   x2 = float(clumpDataBase[i].split()[2])
   y2 = float(clumpDataBase[i].split()[3])
   radius = float(clumpDataBase[i].split()[4])
   eq_radius = float(clumpDataBase[i].split()[5])
   mass = 2670 * (2*z) * (pi * eq_radius**2)
   area += pi * eq_radius**2
   inertia = (pi / coef * 0.315**2 + pi / 2/ coef) * mass * radius**2
   cpos = [(x1+x2)/2.,(y1+y2)/2.,z]
   color = numpy.random.rand(3)
   ID = O.bodies.appendClumped([sphere([x1,y1,z],material=mat,radius=radius,color=color.tolist()), sphere([x2,y2,z],material=mat,radius=radius,color=color.tolist())])
   clump_id.append(ID[0])
   O.bodies[ID[0]].state.pos = cpos
   O.bodies[ID[0]].state.mass = mass
   O.bodies[ID[0]].state.inertia[2] = inertia
   O.bodies[ID[0]].state.blockedDOFs = 'zXYZ'

O.periodic = True
O.cell.setBox(0.35,0.35,z*2)
print len(clump_id)
t_crit = sqrt(2670 * (2*z) * (pi * 0.004**2) / (O.materials[0].young * 0.01578))
O.dt = .5 * t_crit
print O.dt

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

def postConsol():
   for b in O.bodies:
      b.state.blockedDOFs = 'zXY'
      b.state.vel = Vector3.Zero
      b.state.angVel = Vector3.Zero
      b.state.refPos = b.state.pos
      b.state.refOri = b.state.ori
   setContactFriction(0.5)
   O.cell.trsf = Matrix3.Identity
   O.cell.velGrad = Matrix3.Zero
   O.pause()
   biax.doneHook = 'shear()'

def shear():
   print 'Press start to shear..............\n'
   O.pause()
   biax.goal=(-1.e5, -2.e-1, 0)
   biax.stressMask=1
   biax.relStressTol=0.001
   biax.maxUnbalanced=0.05
   biax.maxStrainRate=(0.05, 0.05, 0.)
   biax.doneHook='term()'
   O.engines=O.engines+[PyRunner(command='saveAddData()',iterPeriod=200)]

def saveAddData():
   stress = biax.stress
   strain = biax.strain
   plot.addData(
      strain_x = strain[0],
      strain_y = strain[1],
      stress_x = stress[0],
      stress_y = stress[1]
   )

def term():
   plot.saveDataTxt('biax.txt',vars=('strain_x','strain_y','stress_x','stress_y'))
   O.pause()

O.run()

Why particle along Y axis rotate? How can I fix this problem?
Many thanks for your help.
Best regards, Lifan