Initial stress when I use clump instead of spherical particles

Hello All,

I want to use clumps instead of spherical particles in PeriTriaxController. My problem is at the initial time, the stresses on the wall is more that my goal stress. I don't know where this stresses come from. I will be gratefule if you could help me. Thank you very much.

Best Regards

Alireza

#===================================
from yade import utils, plot
from yade import pack, qt
from datetime import datetime

#==================set a periodic boundary==================

O.periodic=True
O.cell.refSize=(2e-2,2e-2,2e-2)

#================= define the materials =======================

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= True, young=6.81e8, density=1377.5, poisson=0.3, frictionAngle= 0.31, fragile=False, label='aggregate-48'))

#===================clump shape=================================

radz1=[3e-3,3e-3]
poz1= [[-1e-3,0,0],[1e-3,0,0]]
template1= []
template1.append(clumpTemplate(relRadii=radz1,relPositions=poz1))

nums=['t']

mats=['aggregate-48']

coke=(3e-3,4)

color=((0,0,1))

tolerance=[(2e-4)]

temps=[template1]

nums=pack.SpherePack()

nums.makeCloud((0,0,0),(2e-2,2e-2,2e-2),rMean=coke[0],rRelFuzz=tolerance[0],num=coke[1])

O.bodies.append([utils.sphere(c,r,material=mats[0],color=color[0]) for c,r in nums])
O.bodies.replaceByClumps(temps[0],[1.0],discretization=5)

#====================Engine====================
sigmaIso=-1e5

O.engines=[
     ForceResetter(),
     InsertionSortCollider([Bo1_Sphere_Aabb()]),
     InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
     ),
     PeriTriaxController(label='triax',
  goal=(sigmaIso,sigmaIso,sigmaIso),stressMask=7,
  dynCell=True,maxStrainRate=(0.1,0.1,0.1),
  maxUnbalanced=.01,relStressTol=1e-3,
  doneHook='compactionFinished()'
     ),
     NewtonIntegrator(damping=0.4,gravity=[0,0,-10]),

     PyRunner(command='history()',iterPeriod=10,label='recorder'),
]

O.dt=5e-8

def history():
 plot.addData(unbalanced=unbalancedForce(),
  sxx=-triax.stress[0],syy=-triax.stress[1],szz=-triax.stress[2],
  exx=-triax.strain[0],eyy=-triax.strain[1],ezz=-triax.strain[2],
                q=-(triax.stress[2]-triax.stress[1]), p=-(triax.stress[2]+triax.stress[1]+triax.stress[0])/3,
                por=porosity(),i=O.iter,)

        print 'porosity:', porosity(),'Unbalance', unbalancedForce(), 'Strain',-triax.strain[2], 'Stress',-triax.stress[2]