How calculate the repose angle

Asked by Seti on 2015-08-03

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?



Question information

English Edit question
Yade Edit question
No assignee Edit question
Solved by:
Christian Jakob
Last query:
Last reply:
Christian Jakob (jakob-ifgt) said : #1


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:

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

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

O.load(yourSaveFile) #put the name of your save file here

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

height = []
rad_dist = []
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
  rad_dist.append(sqrt( (dist_vec[0])**2 + (dist_vec[1])**2 ))#get radial distance from center axis
  height.append(pos[2]) #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[0][1] - data_for_angle[1][1]) #y entries of picked points (respectively spheres height)
rdiff = abs(data_for_angle[0][0] - data_for_angle[1][0]) #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 ----------

# 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))

Seti (seti) said : #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:
from yade import pack,utils#, qt
pred = pack.inAlignedBox((0,0,0),(20,200,20))
#create material
soil1 = CohFrictMat(young=1e10,poisson=0.2,frictionAngle=radians(0),density=2500.0,normalCohesion=1e6, shearCohesion=80e6,label='soil')
#color=(1,0,0) ----red color
#soil1 = FrictMat(young=1e6,poisson=0.4,frictionAngle=radians(30),density=2500.0,label='soil')
O.materials.append(CohFrictMat(young=1e9,poisson=0.1, frictionAngle = radians(0) , label='wallmat'))
wallmat = O.materials[-1]


             ForceResetter(),#reset forces
                            [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
   # this engine will be called after 20000 steps, only once

# 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
# save the simulation, so that it can be reloaded later, for experimentation
from yade import qt
#from yade import qt
# 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

# 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

Best Christian Jakob (jakob-ifgt) said : #3

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

Seti (seti) said : #4

Thanks Christian Jakob , that solved my question. :)

Seti (seti) said : #5

Thanks Christian Jakob, that solved my question.