Spheres Angle of Repose

Hi all,

I want to calibrate a CohFrictMat in order to achieve an agle of repose of roughly 40°/43° making a heap test with spheres. I tried to tune alphaKr and etaRoll (frictionAngle is fixed at 90°), however I can't see any difference in the result. When the heap is stable (unbalancedForce<0.02) it shows always an angle of about 20°. I ranged both parameters from 1e-03 to 1e2.
Suggestions?

###### MWE ###

# -*- coding: utf-8 -*-

from yade import pack, plot, export, ymport, utils
from yade.gridpfacet import *
from numpy import arctan
import math

RED = '\033[91m'
BLACK = '\033[0m'

#name='_HeapTest'
#vtk_path='VTK'+name
#if not os.path.exists(vtk_path):
# os.makedirs(vtk_path)

O.dt=1e-8
vtk_per=0.1

O.materials.append(CohFrictMat(
                        young=2.85e8,poisson=0.3,density=2600,frictionAngle=radians(90),label='Filling',
                        alphaKr=.2, etaRoll=10, momentRotationlaw=True))
O.materials.append(FrictMat(
                        young=5e8,poisson=0.3,density=7.85e3,frictionAngle=radians(40.),label='FrictMat'))

r3=0.12 #[m]

sp=pack.SpherePack()
sp.makeCloud((.3,-.3,.1),(1.3,.7,2.),rMean=.3*r3,rRelFuzz=0.)
sp.toSimulation()

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),
                       Bo1_GridConnection_Aabb(),
                       Bo1_PFacet_Aabb(),
                       Bo1_Facet_Aabb()]),
     InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),
                        Ig2_Sphere_GridConnection_ScGridCoGeom(),
                        Ig2_Sphere_PFacet_ScGridCoGeom(),
                        Ig2_GridNode_GridNode_GridNodeGeom6D(),
                        Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
                        Ig2_GridConnection_PFacet_ScGeom(),
                        Ig2_PFacet_PFacet_ScGeom(),
                        Ig2_Sphere_Sphere_ScGeom6D(),
                        Ig2_Facet_Sphere_ScGeom6D()],
                      [Ip2_FrictMat_FrictMat_FrictPhys(),
                         Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False)],
                   [Law2_ScGeom_FrictPhys_CundallStrack(),
                         Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True),
                         Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
                         Law2_GridCoGridCoGeom_FrictPhys_CundallStrack()
                        ]),
        GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.9,label='ts'),
        NewtonIntegrator(gravity=(0,0,-9.81),damping=0.8,label='newton'),
        PyRunner(initRun=True,virtPeriod=vtk_per,command='checkUnbalanced()',dead=0),
        PyRunner(initRun=True,virtPeriod=vtk_per,command='addData()'),
        PyRunner(initRun=True,virtPeriod=vtk_per,command='saveData()',label='savedata',dead=1),
]

# BASE FLOOR
kw= {'fixed':True,'material':'FrictMat','color':[1,1,0]}
O.bodies.append(geom.facetBox((1.,0.5,0.),(2.5,2.5,.1),wallMask=1+2+4+8+16,**kw))

def checkUnbalanced():
    if unbalancedForce()<.02:
        export.text('Cumolo.txt')
        O.pause()
        print ''
        print 'TEST ENDED'

def addData():
    plot.addData(iter=O.iter,
                 t=O.time,
                 unbalancedF=unbalancedForce(),
                 )
plot.plots={'t':('unbalancedF')}

O.run()

##### END MWE

Thanks,
Nicola