Using Hertz - Mindlin Contact Law

Hi everybody
I want to carry out triaxial simulation on sand with considering a non-linear contact law (Hertz-Mindlin). I mention the code in below which is working very well with Cohesion moment law but in terms of Hertz- MIndlin there are some problems in the case of creating a pack such as an unacceptable unbalanced (almost 0.67) force for desire void ratio which is around 0.83.
Do you know what should I do to well-create a pack with this non-linear contact law ?
Also, I need to get to know the Hertz-Mindlin equations that are using on YADE, do you know the references of this contact law?

# encoding: utf-8

import matplotlib; matplotlib.rc('axes',grid=True)
import matplotlib.pyplot as plt
from yade import pack
import pylab
from numpy import *
import numpy as np
from math import *

utils.readParamsFromTable(seed=1,num_spheres=10000,compFricDegree = 25.0)
from yade.params import table

seed=table.seed
num_spheres=table.num_spheres# number of spheres
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
confiningS=-1e5
rate=5

psdSizes,psdCumm=[0.00012,0.00014,0.00018,0.00021,],[0.0030,0.039,0.21,1]
psdSizesArray=np.array(psdSizes)
psdSizesMeter=psdSizesArray*0.001 #Convert the size of particles to meter
sp=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(0.0001,0.0001,0.0001)
sp.makeCloud(minCorner=mn,maxCorner=mx,psdSizes=psdSizesMeter,psdCumm=psdCumm,distributeMass=True,num=num_spheres,seed=seed)
sp.psd(bins=10,mass=True)

## create material #0, which will be used as default
O.materials.append(FrictMat(young=5e9,poisson=.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=5e9,poisson=.3,frictionAngle=0,density=0,label='frictionless'))

## create walls around the packing
walls=aabbWalls((mn,mx),thickness=0,material='frictionless')
wallIds=O.bodies.append(walls)

O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

triax=TriaxialStressController(
 internalCompaction=False,
 goal1=confiningS,
 goal2=confiningS,
 goal3=confiningS,
 max_vel=10,
 label="triax"
)

newton=NewtonIntegrator(damping=0)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_MindlinPhys()],
  [Law2_ScGeom_MindlinPhys_Mindlin(includeAdhesion=False,includeMoment=False)]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 VTKRecorder(iterPeriod=1000,fileName="./export",recorders=['spheres','intr','color'],dead=1,label='VTK'),
 triax,
 newton
]

#s11=[];e22=[]
while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  print ('unb=', unb,'meanstress:',abs(triax.goal1-triax.meanStress)/abs(triax.goal1),'e=', triax.porosity/(1-triax.porosity), 'Size=', (triax.width))
  if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.01 and (triax.porosity)/(1-triax.porosity)<0.8399999999 :
    break
print ("Particle.Distribution.is.stable")
 ########################################################
 ##Turn into triaxil test
 ########################################################

finalFricDegree = 18.0
triax.internalCompaction=False
setContactFriction(radians(finalFricDegree))
while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.01:
    break
    #ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2]
    #si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2]

triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width
triax.stressMask=5
triax.goal1=triax.goal3=confiningS
triax.goal2=-rate
triax.wall_bottom_activated=0
#O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True
VTK.dead=1
#cohesiveLaw.always_use_moment_law=False
#recorder.dead=0

Exaxil=[5.98e-04,6.11e-04,0.00125,0.0013,0.00249,0.0031,0.00491,0.00729,0.010,0.017,0.025]
Exdev=[8.69,16.30,51.09,81.52,92.40,110.88,142.40,170.67,196.77,218.54,229.44]

ee=[];e22=[];s22=[]
file=open('Results'+'.txt',"w")

while 1:
  O.run(1000,True)
  ee.append((triax.porosity)/(1-triax.porosity))
  s22.append(((-triax.stress(triax.wall_top_id)[1])-(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0)*1e-3)
  e22.append((-triax.strain[1]))
  print('ax=',(-triax.strain[1]), 'deviatoric=',(((-triax.stress(triax.wall_top_id)[1])-(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0)*1e-3), 'ee=', (triax.porosity)/(1-triax.porosity), 's11=',((-triax.stress(triax.wall_top_id)[1])*1e-3), 's33=',(((-triax.stress(triax.wall_right_id)[0])-(triax.stress(triax.wall_front_id)[2]))/2.0)*1e-3)
  if abs(-triax.strain[1])>0.1 :
   break
print(" Triaxil test is finished ")

file.write(str(s22)+" "+str(e22)+"\n")
file.write(str(s22)+" "+str(ee)+"\n")
file.close()

plt.figure()
plt.xlabel('Axial strain')
plt.ylabel('deviatoric')
plt.grid(True)
plt.plot(Exaxil,Exdev,linestyle='-',marker="*", markersize=8,color='g',linewidth=2)
plt.plot(e22,s22,linestyle='-',marker="^", markersize=8,color='r',linewidth=2)
plt.legend(('Exprimental data','simulation data'),loc='upper right', shadow=True)

plt.show()