Porosity and friction angle

Hello
Yade version: 1.20.0
I want to simulate a cylinder filled with spherical particles . After reaching to specific unbal force cylindrical confining wall is going to move upward with specific velocity. I am going to look at the height of the particles after several steps with different material properties.
My question: I am trying to change porosity of materials by changing friction angle and reach to a specific one. However, changing friction angle does not effect porosity in my model. I get very close porosity value by changing the friction angle. (even sometimes by decreasing the friction angle, porosity increases). Also porosity I am getting is not close to my scope( what I am getting in PFC). I am wondering what can be changed in my code for solving it?
 I am posting my code below. Also I used [1] for porosity function in my code as you advised others with similar problem.

[1]:https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py#L141
#######################################################################################
# import yade modules that we will use below
from yade import pack,utils, plot,ymport,qt
#######################################################################################

fr1 = 45;fr2 = 45;rho=2650.0

Mat1 = O.materials.append(ViscElMat(kn=10.0e4,ks=10.0e4,frictionAngle=radians(fr1),density=rho))
Mat2 = O.materials.append(ViscElMat(kn=10.0e4,ks=10.0e4,cn=0.1,cs=0.1, frictionAngle=radians(fr2),density=rho))

#######################################################################################

# create Cylinder from facets
Cylinder_mo=O.bodies.append(geom.facetCylinder((0,0,5),radius=1.0,height=10.0,segmentsNumber=15,wallMask=4,material=Mat1))
Wall_below=O.bodies.append(utils.wall(position=(0,0,0),axis=2,sense=0,material=Mat1))
# create empty sphere packing
# sphere packing is not equivalent to particles in simulation, it contains only the pure geometry
sp=pack.SpherePack()
# generate randomly spheres with uniform radius distribution
sp.makeCloud((-0.5,-0.5,0.0),(0.5,0.5,4.0),rMean=.05)#,material=Mat1)
sphe_particles=O.bodies.append([sphere(c,r,material=Mat2) for c,r in sp])

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      # handle sphere+sphere and facet+sphere collisions
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),

   PyRunner(command='hightcontrol()',iterPeriod=150000),

    NewtonIntegrator(gravity=(0,0,-9.81),damping=0.0),
    # call the addPlotData function every 200 steps
     PyRunner(command='addPlotData()',iterPeriod=100),
   ]

O.dt=.5*PWaveTimeStep()

# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.energy['energyName'] subsequently
O.trackEnergy=True

# if the unbalanced forces goes below .05, the packing
# is considered stabilized, therefore we stop collected
# data history and stop
####################################################################
#def checkUnbalanced():
  # if unbalancedForce()<.01:
     # O.pause()
     # plot.saveDataTxt('bbb.txt.bz2')
      #plot.saveGnuplot('bbb') is also possible
####################################################################
####################################################################
def hightcontrol():
 for b in sphe_particles:
   if O.bodies[b].state.pos[2]>0.20:
      O.bodies.erase(b);

def Rot_resist():
 for b in sphe_particles:
    O.bodies[b].state.blockedDOFs='XYZ'
    O.bodies[b].state.angVel=0.99*O.bodies[b].state.angVel

####################################################################

# collect history of data which will be plotted
def addPlotData():
   # each item is given a names, by which it can be the unsed in plot.plots
   # the **O.energy converts dictionary-like O.energy to plot.addData arguments
   plot.addData(i=O.iter,unbalanced=unbalancedForce())#,**O.energy)

plot.plots={'i':('unbalanced')}#,None,O.energy.keys)}

# show the plot on the screen, and update while the simulation runs
plot.plot()
O.run(145000);O.wait() # wait at this step then go to next line .
O.save('Mymodel2.yade')

#porosity check by changing friction angle
import sys
while utils.porosity()>0.45:
    fr2=0.95*fr2
    setContactFriction(radians(fr2))
#vo=utils.aabbDim(cutoff=0.0,centers=False)
    print "porosity:",utils.porosity()," Fricition:",fr2
    #print "porosity:",volume
    sys.stdout.flush()
    O.run(5000);O.wait()
    O.save('Mymodel3.yade')

O.dt=.5*PWaveTimeStep()