Shear loading

I made a simple model to understand what’s going on. I have two spheres:
sphere(center=(0,0,0),radius=0.5,fixed=True,material='mat1'),
sphere((0,0,1),radius=0.5,material='mat1'),
The following is the boundary condition:
O.bodies[1].state.blockedDOFs = "xyzXYZ"
O.bodies[1].state.vel = (0.0001,0,0)
O.bodies[1].state.angVel = (0,0,0)

I am applying a velocity in the x direction on sphere #1 which should cause a shear reaction force to occur. The following are my questions:
1. When I track the forces on sphere #1, I notice that there is a force in the x-direction (which is expected) and a force in the z-direction (uniaxial). Where is the uniaxial force coming from if I only applied a shear boundary condition? The contrary example to this would be to apply a uniaxial loading condition which would result in only a uniaxial force.
2. I have read [1] and [2] thoroughly. I understand that the tangential plane rotates as the spheres shear and that the shear force is calculated incrementally. Is deltaUs in [2] the change in tangential displacement with respect to time or is it from the initial configuration?
3. I tried to do a simple DEM model on MATLAB where I had two 2D spheres and I applied the same loading condition as my YADE model while calculating the uniaxial and shear forces on sphere#1. I could not get the same results as YADE predicts. Is there a document that explains clearly the steps. I read [1], but I want something with more detail.
4. My code is running well, so I assume that the my code is correct.

[1] https://yade-dem.org/doc/formulation.html#fig-shear-2d
[2] A DEM model for soft and hard rocks: Role of grain interlocking on strength

Thank you soo much for your time.

MY CODE:
# MATERIAL PROPERTIES
intR=20 # allows near neighbour interaction (can be adjusted for every packing)
DENS=2500 # Density
YOUNG=1800
FRICT=7
ALPHA=0.1
TENS=100
COH=100
iterMax = 100

O.materials.append(JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH,label='mat1'))

# PARTICLES
O.bodies.append([
 # 0
 sphere(center=(0,0,0),radius=0.5,fixed=True,material='mat1'),
 # 1
 sphere((0,0,1),radius=0.5,material='mat1')
])

nbSpheres = 2;

# BOUNDARY CONDITION
O.bodies[1].state.blockedDOFs = "xyzXYZ"
O.bodies[1].state.vel = (0.0001,0,0)
O.bodies[1].state.angVel = (0,0,0)
# FUNCTIONAL COMPONENTS
# 1
def addPlotData1():
        plot.addData(t=O.time,DX1 = O.bodies[1].state.pos[0],DY1 = O.bodies[1].state.pos[1],DZ1 = O.bodies[1].state.pos[2], FX1 = O.forces.f(1)[0], FY1 = O.forces.f(1)[1], FZ1 = O.forces.f(1)[2])

# 2
def addPlotData2():
 plot.saveDataTxt('/home/nabid/Desktop/SCRIPTS/dataFile.txt',vars=('t','DX1','DY1','DZ1','FX1','FY1','FZ1'))

#PRINTING
from yade import plot
from pprint import pprint

# SIMULATION LOOP
O.engines=[
 PyRunner(command='addPlotData1()',iterPeriod=1),
 PyRunner(command='addPlotData2()',iterPeriod=1),
     ForceResetter(),
     InsertionSortCollider([Bo1_Sphere_Aabb()]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR)],
    [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
    [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
  ),
     GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5, defaultDt=utils.PWaveTimeStep()),
 NewtonIntegrator(damping=0.1),
]

#TIME STEP
O.dt=0.5e-4*PWaveTimeStep()

#PLOTTING
plot.plots={'DX1':('FX1'),'DX1':('FZ1')}
plot.plot()

# SAVE SIMULATION