# How calculate the repose angle

I just have written a simple script which just models the falling down of particles because of their weight.. The end of simulation, I need to calculate the repose angle of the accumulated grains ( the angle of triangle which is shaped because of falling grains due to gravity). but I do not have any idea in this issue, would you please instruct me?

Cheers

Seti

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Christian Jakob
Solved:
2015-08-10
Last query:
2015-08-10
2015-08-06
 Christian Jakob (jakob-ifgt) said on 2015-08-04: #1

Hi,

I had the same problem some years ago. I used pythons scatter function for this purpose.
You need to give some input data at the beginning of the script and it should work (not tested).
A scatter-plot window should open and you can get the angle by clicking at two points in the window.
Here is the script:

#!/usr/bin/python
# -*- coding: utf-8 -*-

center_model = Vector3(0,0,0) #put the center of your model here

outputFilename = 'yourOutputName.out' #put your output filename here

height = []
for ii in range(0,len(O.bodies)):
b = O.bodies[ii]
if b and isinstance(b.shape,Sphere):
pos = b.state.pos
dist_vec = pos - center_model
height.append(pos) #get height and write into list

from pylab import *

data_for_angle = ginput(0,0) #get the coordinates by clicking two points, then click on the mouse wheel

zdiff = abs(data_for_angle - data_for_angle) #y entries of picked points (respectively spheres height)
rdiff = abs(data_for_angle - data_for_angle) #x entries of picked points (respectively spheres radial distance from center)

alpha = atan(zdiff/rdiff)
alpha = 180*alpha/math.pi #in [°]

'''
z |-
d | -
i | -
f | ( - <- tan(alpha) = zdiff / rdiff
f ----------
rdiff
'''

# print result and write output:
print 'angle of repose is ',alpha,' degree'

f = open(outputFilename,'a')
f.write('result for %s: %.2f degree\n' % (part,alpha))
f.close()
exit()

 Seti (seti) said on 2015-08-06: #2

Hi Christian,

Thank you for your reply, I have run your script. As you have mentioned a scatter-plot window opens however by clicking at points nothing happen.Would you please instruct me in this issue?

Here is my script:
#!/usr/bin/python
pred = pack.inAlignedBox((0,0,0),(20,200,20))
#create material
#color=(1,0,0) ----red color
O.materials.append(soil1)
O.bodies.append(utils.wall(0,axis=1,sense=1))
O.materials.append(CohFrictMat(young=1e9,poisson=0.1, frictionAngle = radians(0) , label='wallmat'))
wallmat = O.materials[-1]

spheres=SpherePack()
spheres.toSimulation()
#O.bodies.append(spheres)
#

#
O.engines=[
ForceResetter(),#reset forces
InsertionSortCollider([Bo1_Wall_Aabb(),Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Wall_Sphere_ScGeom()], # collision geometry
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()], # collision "physics"
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()] # contact law -- apply forces
),
# apply gravity force to particles
# damping: numerical dissipation of energy
NewtonIntegrator(damping=0.5,gravity=(0,-9.81,0)),
#qt.SnapshotEngine(fileBase='3d-',iterPeriod=200,label='snapshot'),
# this engine will be called after 20000 steps, only once
#PyRunner(command='finish()',iterPeriod=20000)
]

# set timestep to a fraction of the critical timestep
# the fraction is very small, so that the simulation is not too fast
# and the motion can be observed
O.dt=1*utils.PWaveTimeStep()
#makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=10000)
# save the simulation, so that it can be reloaded later, for experimentation
#O.saveTmp()
O.save('Modified.txt.bz2')
qt.View()
#O.run()
#qt.View()
#O.run()
# this function is called when the simulation is finished
#def finish():
# snapshot is label of qt.SnapshotEngine
# the 'snapshots' attribute contains list of all saved files
#makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=10000)
#O.pause()

# set parameters of the renderer, to show network chains rather than particles
# these settings are accessible from the Controller window, on the second tab ("Display") as well
#rr.shape=False Christian Jakob (jakob-ifgt) said on 2015-08-06: #3

I think you need to click on the mouse wheel ... (as mentioned in the script)

 Seti (seti) said on 2015-08-10: #4

Thanks Christian Jakob , that solved my question. :)

 Seti (seti) said on 2015-08-10: #5

Thanks Christian Jakob, that solved my question.