Difference between CundallStrack and HertzMindlin

Hi,
I'm doing uniaxial compression of particles in cylinder container. I'm comparing 2 contact models (Hertz Mindlin and Cundall Strack). I noticed that the load needed to obtain volume ratio (volume of particles/volume of container) 1 is much higher when using Cundall Strack model (12 000 N, comparing to 2000 N when using Hertz Mindlin). All the parameters are the same in both simulations. Does anyone know why this happens? Also the load in both cases increases linearly, and by experiments after some volume ratio it should increase kind of exponentionally (for example, the steep of the load curve between volume ratios 0.85 and 1 is way steeper than for less volume ratios).
Thanks in advance

MWE:

from yade import pack, plot, export
import math
#sphere parameters
r=0.0025
rr=0.05
n=800

#material definition
E=31e6
v=0.3897
Ro=1140
dmp=0.05
fr=15
O.materials.append(FrictMat(young=E, poisson=v, frictionAngle=radians(fr), density=Ro, label='spheres'))

Ec=2.1e9
vc=0.33
Roc=7850
frc=22
O.materials.append(FrictMat(young=Ec, poisson=vc, frictionAngle=radians(frc), density=Roc, label='walls'))

#cylinder container
R=0.015
h=10*R
c=(0,0,h/2)
cylinder = geom.facetCylinder(center=c, radius=R, height=h, segmentsNumber=12, wallMask=6, material='walls')
O.bodies.append(cylinder)

#funnel definition
c1=(0,0,h)
dB=8*R
dO=2*R
hB=1.5*h
hO=0.5*h
hP=0
funnel = geom.facetBunker(center=c1, dBunker=dB, dOutput=dO, hBunker=hB,hOutput=hO, hPipe=hP, segmentsNumber=20, wallMask=4, material='walls')
O.bodies.append(funnel)

#sphere packing
sp = pack.SpherePack()
sp.makeCloud((-dB/3,-dB/3,1.5*h),(dB/3,dB/3,1.9*h), rMean = r, rRelFuzz = rr, num = n)
sp.toSimulation(material = 'spheres')

#setting simulation engines
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb(), Bo1_Wall_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom(), Ig2_Wall_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_MindlinPhys()],
                [Law2_ScGeom_MindlinPhys_Mindlin(label='Mindlin')]
        ),
        NewtonIntegrator(gravity=(0, 0, -9.81), damping=dmp),
        PyRunner(command='checkUnbalanced()', realPeriod=2, label='checker'),
]

#timestep
O.dt = 0.3*PWaveTimeStep()

#checking the status of gravity deposition and loading if done
def checkUnbalanced():
 if O.iter < 150000:
  return
 if unbalancedForce()>1:
  return
 #adding the loading plate
 O.bodies.append(wall(max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)]), axis=2, sense=-1))
 global plate
 plate = O.bodies[-1]
 w=plate.state.pos[2]
 plate.state.vel =(0,0,-0.05)
 O.engines = O.engines + [PyRunner(command='addPlotData()', iterPeriod=100)]
 checker.command = 'unloadPlate()'

#checking the loading status and unloading if done
def unloadPlate():
 V1=((2*R)*(2*R)*3.14/4)*plate.state.pos[2]
 V2=n*4*3.14*r*r*r/3
 Vr=V2/V1
 Fz = O.forces.f(plate.id)[2]
 S=Fz/(4*R*R*3.14/4)
 if Vr>=1.1:
  O.pause()

#plotting the data
def addPlotData():
 if not isinstance(O.bodies[-1].shape, Wall):
  plot.addData()
  return
 Fz = O.forces.f(plate.id)[2]
 w=plate.state.pos[2]
 V1=((2*R)*(2*R)*3.14/4)*w
 V2=n*4*3.14*r*r*r/3
 Vr=V2/V1
 S=Fz/(4*r*r*3.14/4)
 plot.addData(Fz=Fz, w=plate.state.pos[2], unbalanced=unbalancedForce(), Vr=V2/V1, S=Fz/(4*R*R*3.14/4))

plot.plots = {'w':('Fz'),'Vr':('S'),}
plot.plot()

O.run()