hertz model

Hi All,

I'm trying to use this constitutive model (Law2_ScGeom_MindlinPhys_Mindlin()) to simulate the triaxial tests.
here is my code:
###################################
############################
### DEFINING PARAMETERS ###
############################
from yade import pack, plot
targetPorosity=0.3682
compFricDegree=45.0
finalFricDegree=19.5
rate=-0.5
damp=0.7
stabilityThreshold=0.001
young=4e8
mn,mx=Vector3(0,0,0),Vector3(4e-3,8e-3,4e-3)
##################################################################
###################### material for sphere and wall ##############
##################################################################
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2648,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
psdSizes=[0.12e-3,0.14e-3,0.16e-3,0.18e-3,0.19e-3,0.20e-3,0.22e-3,0.24e-3,0.29e-3,0.38e-3]
psdCumm=[0,11.1,22.2,33.3,44.4,55.5,66.6,77.7,88.8,100]
sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,num=11526,distributeMass=1,seed=1)
sp.toSimulation(material='spheres')
O.dt=.5*PWaveTimeStep()
############################
### DEFINING ENGINES ###
############################
triax=TriaxialStressController(
 maxMultiplier=1.0+2e5/young,
 finalMaxMultiplier=1.0+2e4/young,
 thickness=0,
 stressMask=7,
 internalCompaction=False
)
newton=NewtonIntegrator(damping=damp)
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
   [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_FrictMat_FrictMat_MindlinPhys( betan =0.2 , betas = 0.2, krot = 0.2 , eta = 0.5 )],
   [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom_MindlinPhys_Mindlin()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8),
 triax,
 #VTKRecorder(fileName='george-triaxial-',recorders=['all'],iterPeriod=1000),
 #TriaxialStateRecorder(iterPeriod=2000,file='WallStresses-george'+table.key),
 newton
]
########################################
#### APPLYING CONFINING PRESSURE ###
########################################
triax.goal1=triax.goal2=triax.goal3=-100000
while 1:
 O.run(3000,True)
 unb=unbalancedForce()
 print 'unbalanced force: ',unb,'mean stress: ',triax.meanStress
 if unb<stabilityThreshold and abs(-100000-triax.meanStress)/-100000<0.001:
  break
####################################################
#### REACHING A SPECIFIED POROSITY PRECISELY ###
####################################################
##import sys #this is only for the flush() below
while triax.porosity>targetPorosity:
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 O.run(500,1)
#set new friction angle
triax.wall_bottom_activated=True
triax.wall_top_activated=True
triax.wall_left_activated=True
triax.wall_right_activated=True
triax.wall_back_activated=True
triax.wall_front_activated=True
setContactFriction(radians(finalFricDegree))
while 1:
 O.run(3000,True)
 unb=unbalancedForce()
 print 'unbalanced force: ',unb,'mean stress: ',triax.meanStress
 if unb<stabilityThreshold and abs(-100000-triax.meanStress)/-100000<0.001:
  break
print "porosity= ",triax.porosity
triax.goal1=triax.goal2=triax.goal3=-100000
while 1:
    O.run(1000, True)
    #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
    unb=unbalancedForce()
    print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
    if unb<stabilityThreshold and abs(-100000-triax.meanStress)/-100000<0.001:
        break
#O.save('compactedState'+key+'.yade.gz')
print "### Compacted state saved ###"
print "porosity= ",triax.porosity
re11=-triax.strain[0]
re22=-triax.strain[1]
re33=-triax.strain[2]
################################
##### DEVIATORIC LOADING ###
################################
triax.stressMask = 5
newton=NewtonIntegrator(damping=0.0)
#now goal2 is the target strain rate
triax.goal2=rate
# we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-100000
triax.goal3=-100000
triax.strainRate2=1
triax.strainRate1=triax.strainRate3=1000.0
######################################################
#### Example of how to record and plot data ###
######################################################
def history():
   plot.addData(e11=-triax.strain[0]-re11, e22=-triax.strain[1]-re22, e33=-triax.strain[2]-re33,
        ev=-triax.strain[0]-re11-triax.strain[1]-re22-triax.strain[2]-re33,
      s11=-triax.stress(triax.wall_right_id)[0],
      s22=-triax.stress(triax.wall_top_id)[1],
      s33=-triax.stress(triax.wall_front_id)[2],
             i=O.iter)
 plot.saveDataTxt('calibration-george.txt')
if 1:
  ## include a periodic engine calling that function in the simulation loop
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=1000,command='history()',label='recorder')]+O.engines[5:7]
  #O.engines=O.engines[0:5]+[PyRunner(iterPeriod=2000,command='savedata()',label='recorder_1')]+O.engines[5:7]
  ##O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder'))
else:
  O.engines[4]=PyRunner(iterPeriod=1000,command='history()',label='recorder')
  #O.engines[4]=PyRunner(iterPeriod=2000,command='savedata()',label='recorder_1')
O.run(1000000,True)
####################################################
I'm using High-performance computing with 22 CPUs.
########################################
There are no errors.
#######################
The total number is around 10000.
The unbalanced force did not reach the threshold (= 0.001) after one day.
I want to know if this situation is reasonable? (because when I use this constitutive model: Law2_ScGeom6D_CohFrictPhys_CohesionMoment, the unbalanced force reaches the threshold just 7 hours ).
#############################

best,
Yong