Triaxial cohesive

Dear all,

I have done some small changes on cohesive triaxial script - changing the normal and cohesion values for two type of materials - but the results does not make sense - unbalace force nan and mean stress zero

can you please advise me what is my mistake thanks

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division

from yade import plot,pack,timing
import time, sys, os, copy
from yade import ymport, utils, pack, export
#from pylab import *
import math
import numpy as np

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################
# default parameters or from table
readParamsFromTable(noTableOk=True, # unknownOk=True,
 targetPorosity=0.45, #the porosity we want for the packing
 num_spheres=3000,# number of spheres

 young=236e9,
 poisson=.2,
 density=2600,
 frictionAngle=radians(30),
 normalCohesion=3.e6,
 shearCohesion=1.e6,
 etaRoll=0.1,

 compFricDegree=30, # contact friction during the confining phase
 finalFricDegree=30, # contact friction during the deviatoric loading
 key='_triax_base_', # put you simulation's name here

 rate=-0.02, # loading rate (strain rate)
 damp=0.7, # damping coefficient
 stabilityThreshold=0.01, # we test unbalancedForce against this value in different loops (see below)
)
from yade.params import table
from yade.params.table import *
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing
## create materials for spheres and plates
mat1=CohFrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,normalCohesion=normalCohesion,shearCohesion=shearCohesion,momentRotationLaw=True,etaRoll=etaRoll,label='cement')
mat_1=O.materials.append(mat1)
mat2=CohFrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,normalCohesion=3.e6,shearCohesion=1.e6,momentRotationLaw=True,etaRoll=etaRoll,label='glass_sphere')
mat_2=O.materials.append(mat2)
O.materials.append(FrictMat(young=2*young,poisson=.25,frictionAngle=0,density=0,label='frictionlessWalls'))

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

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
O.bodies.append([sphere(center,rad,material='cement') for center,rad in sp])
cement=[]
glass_sphere=[]
for b in O.bodies:
 if not isinstance(b.shape,Sphere): # change material only on spheres
    continue
 if random.random() < 0.5:
  b.mat = mat1
  b.shape.color = (1,0,0)
  b.state.mass*=mat1.density/mat1.density
  cement.append(b.id)
 else:
  b.mat = mat2
  b.shape.color = (0,1,1)
  b.state.mass*=mat2.density/mat1.density
  glass_sphere.append(b.id)
############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
 maxMultiplier=1.+2.4e5/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2.4e4/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask = 7,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  #box-sphere interactions will be the simple normal-shear law, we use ScGeom for them
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
  #Boxes will be frictional (FrictMat), so the sphere-box physics is FrictMat vs. CohFrictMat, the Ip type will be found via the inheritance tree (CohFrictMat is a FrictMat) and will result in FrictPhys interaction physics
  #and will result in a FrictPhys
  [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")],
  #Finally, two different contact laws for sphere-box and sphere-sphere
  [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
   useIncrementalForm=True, #useIncrementalForm is turned on as we want plasticity on the contact moments
   always_use_moment_law=False, #if we want "rolling" friction even if the contact is not cohesive (or cohesion is broken), we will have to turn this true somewhere
   label='cohesiveLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]
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
print 'timestep at confinedState:',O.iter

import sys #this is only for the flush() below
while triax.porosity>targetPorosity:
 ## we decrease friction value and apply it to all the bodies and contacts
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 #print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 ## while we run steps, triax will tend to grow particles as the packing
 ## keeps shrinking as a consequence of decreasing friction. Consequently
 ## porosity will decrease
 O.run(1000,1)
print 'porosity:',triax.porosity
print 'timestep at compactedState:',O.iter
print "### Compacted state saved ###"

##############################
### 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))

##set stress control on x and z, we will impose strain rate on y
triax.stressMask = 5
##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

newton.damping=0.1
O.engines[2].lawDispatcher.functors[1].always_use_moment_law = True
#We assign cohesion to all contacts at the next iteration
O.engines[2].physDispatcher.functors[1].setCohesionNow = True

O.run(1000,True)
O.saveTmp()

#####################################################
### Example of how to record and plot data ###
#####################################################

#from yade import plot

### a function saving variables
def history():
   plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
        ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
      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)

if 1:
  ## include a periodic engine calling that function in the simulation loop
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
  ##O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder'))
else:
  ## With the line above, we are recording some variables twice. We could in fact replace the previous
  ## TriaxialRecorder
  ## by our periodic engine. Uncomment the following line:
  O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder')

O.run(1000,True)

### declare what is to plot. "None" is for separating y and y2 axis
#plot.plots={'i':('e11','e22','e33',None,'s11','s22','s33')}
### the traditional triaxial curves would be more like this:
#plot.plots={'e22':('s11','s22','s33',None,'ev')}
plot.plots={'e22':('s11','s22')}

## display on the screen (doesn't work on VMware image it seems)
plot.plot()

##### PLAY THE SIMULATION HERE WITH "PLAY" BUTTON OR WITH THE COMMAND O.run(N) #####

## In that case we can still save the data to a text file at the the end of the simulation, with:
plot.saveDataTxt('400,alphaKr=0.6,normalCohesion=1e9,shearCohesion=.3e9,young=236e9,j'+key)
##or even generate a script for gnuplot. Open another terminal and type "gnuplot plotScriptKEY.gnuplot:
plot.saveGnuplot('400,alphaKr=0.6,normalCohesion=1e9,shearCohesion=.3e9,young=236e9,j'+key)
rr=yade.qt.Renderer()
rr.shape=True
rr.intrPhys=True

Thanks