Error in activating Moment (includeMoment=True) in Mindlin contact law

Dear all,

I am trying to model a drained triaxial test in dry sand. I used Mindlin contact law for contact between spheres. Also, I want to consider bending moment and activate includeMoment=True in Law2_ScGeom_MindlinPhys_Mindlin(). But i couldn't make a pack or continue shear loading after making pack and consolidation(I see Error nan in terminal). I need to model rolling resistance to get experimental shear strength. Therefore, i gave some value to eta for considering plastic bending moment.

Would you please look at my script or run it and inform me about mistakes if there are?

Thank you for your help,
Elham

# 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=1000,compFricDegree =3.0)
from yade.params import table

seed=table.seed
num_spheres=table.num_spheres
compFricDegree = table.compFricDegree
confiningS=-50000
rate=1

Young=200 #MPa
Poisson=0.3

## creat a packing with a specific particle side distribution (PSD)
psdSizes,psdCumm=[0.075,0.107,0.153,0.192,0.211,0.232,0.252,0.279,0.303,0.462],[0.035,0.065,0.246,0.5,0.6,0.7,0.8,0.9,0.965,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.002,0.004,0.002)
sp.makeCloud(minCorner=mn,maxCorner=mx,psdSizes=psdSizesMeter,psdCumm=psdCumm,distributeMass=True,num=num_spheres,seed=seed)

O.materials.append(FrictMat(young=20e7,poisson=.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=20e7,poisson=.3,frictionAngle=0,density=0,label='frictionless'))

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=True,
 goal1=confiningS,
 goal2=confiningS,
 goal3=confiningS,
 max_vel=10,
 label="triax"
)

newton=NewtonIntegrator(damping=0.4)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_MindlinPhys(en=0.9,es=0.9,eta=0.8,krot=1)],
  [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True,label='mindlin')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]

while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  print ('Porosity:',triax.porosity,'unb:',unb,'meanstress:',abs(triax.goal1-triax.meanStress)/abs(triax.goal1))
  if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
    break

print ("Particle Distribution is Stable")

print ('Porosity:',triax.porosity)

#############################
## REACH NEW EQU. STATE ###
#############################
finalFricDegree = 26.5 # contact friction during the deviatoric loading

#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False

# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  print ('Porosity:',triax.porosity,'unb:',unb,'meanstress:',abs(triax.goal1-triax.meanStress)/abs(triax.goal1))
  if unb<0.001 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
    break

triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width

print ("Time for the deviator loads")

print ('Porosity:',triax.porosity)

##########################
## Start test ###
##########################

AxialstrainMatrix1=[]
DeviatoricStressMatrix1=[]
VolumetricStrainMatrix1=[]
PrincipalStressRatioMatrix1=[]
ESecanti=[]

triax.stressMask=5
triax.goal2=-rate
triax.goal1=confiningS
triax.goal3=confiningS

key='E_'+str(Young)+'_noo_'+str(Poisson)

file=open('Results'+key+'.txt',"w")
while 1:
 O.run(1000,True)
 AxialstrainMatrix1.append((-triax.strain[1]))
 DeviatoricStressMatrix1.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)
 PrincipalStressRatioMatrix1.append(abs((-triax.stress(3)[1])/((-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0)))
 VolumetricStrainMatrix1.append((triax.strain[0]+triax.strain[1]+triax.strain[2]))
 ESecanti.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/(-triax.strain[1]))
 print ('eyy=',triax.strain[1],'Sxx=', triax.stress(0)[0],'Syy=',triax.stress(3)[1], 'Szz=',triax.stress(4)[2],'e=', (triax.porosity)/(1-triax.porosity))

 if triax.strain[1] <-0.03:
  break
file.write(str(AxialstrainMatrix1)+" "+str(DeviatoricStressMatrix1)+" "+str(PrincipalStressRatioMatrix1)+" "+str(VolumetricStrainMatrix1)+" "+"\n")
file.close()