Results of Cohesive triaxial - match maker

Hi all,

Back to my previous question #663085 - I have updated the script as per below - it is working without any error however when I change the percentage of materials from 10% to 90%( through this line :- if random.random() < 0.9:) the stress- strain results does not change at all.
Even I am not sure if "DEVIATORIC LOADING" phase is impaling, because the stress value does not go more than 100kpa in all different ratios of the materials combinations.

Would you please advise me where is the issue?
Thanks
Seti

####
#!/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.e3,shearCohesion=1.e3,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.9:
  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"), Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(normalCohesion=MatchMaker(matches=((mat_1,mat_1,3.e6),(mat_1,mat_2,3.e4),(mat_2,mat_2,3.e3))),shearCohesion=MatchMaker(matches=((mat_1,mat_1,1.e6),(mat_1,mat_2,1.e4),(mat_2,mat_2,1.e3))))],
  #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),
 TriaxialStateRecorder(iterPeriod=1000,file='90%,matchmaker,j'+key),
 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)
##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
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(100,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('90%,matchmaker '+key)
##or even generate a script for gnuplot. Open another terminal and type "gnuplot plotScriptKEY.gnuplot:
plot.saveGnuplot('90%,matchmaker'+key)
rr=yade.qt.Renderer()
rr.shape=True
rr.intrPhys=True