Hello,everyone.
I would like to ask if the dynamic Angle of repose can be output at the moment of rotation of the drum. Need to add code?
Hi Xutong. This is another case where youāre going to have to use some UDF (userdefined functions) to compute the quantity you are interested in. If you come up with a good implementation for AoR, you can share it here.
Searching the MFiX forum, there have been a few mentions of Angle of Repose in the past, but thereās no working examples.
Good luck!
ā Charles
On further reflection, you might be able to use a monitor for this: perhaps monitoring the mean particle location (i.e. center of mass of particles) would suffice.
Just look at how far from the vertical the center of mass is. A simple arctan of the (centerofparticlemass  centerofcylinder) should get you the AoR
āCharles
Unfortunately, there is no monitor for particle position. youāre going to have to write a UDF.
I made a very simple implementation of this idea in the rotating_drum tutorial. Hereās the changes I made to usr1_des.f
:
 /home/cgw/Work/NETL/mfix/tutorials/dem/rotating_drum_3d/usr1_des.f 20210728 07:45:00.972978590 0500
+++ ./usr1_des.f 20211116 13:42:54.119068198 0600
@@ 35,11 +35,17 @@
use discretelement
use run
use usr

+ use constant, only: pi
+
implicit none

+
+ double precision :: mean(3)
integer :: i
+ mean = sum(des_pos_new, 1) / max_pip
+ write(*,*) "MEAN POS", mean
+ write(*,*) "ANGLE", 180*atan2(mean(2), mean(1))/PI + 90
+
! Interpolate data only in root process
if(mype == pe_io) then
call interpolate_keyframe_data(time)
The variables max_pip
and des_pos_new
are defined in the module discreteelement
, for which we already have a use
statement, so they are already accessible.
All this does is compute the center of mass of all the particles, and then how far from the vertical axis that is. When the particles are settled, this should match the angle of repose. But note that this computation will always return some angle, even when the angle of repose is not welldefined (like at the beginning of the simulation before all the particles have dropped). Itās up to you to decide if this algorithm is good enough or if you need some fancier way to compute angle of repose. This code is just to test the idea  it doesnāt handle DMP, for example, and the data needs to be saved in a better format (this just writes to the console). Hereās a plot of the reported angle during the run of rotating_drum
(the drum reverses direction part way through)
Please let me know if this is useful,
ā Charles
Hereās another run, with the direction reversal removed. The drum rotates at a constant rate for the first 8 seconds of the simulation, then is stopped.
You can see that the computed AoR oscillates for a bit at the beginnng, then takes quite a while to settle down. This could be because the algorithm is sensitive to statistical outliers, or perhaps itās a real effect  it may take a few seconds for the balls to find their optimum packing and to lose their kinetic energy from the initial drop.
Youāll have to decide if this centerofmassbased approach is good enough, or if you have to do some more complex measurements to determine AoR.
@Zhangxutong 
hereās a much better result. In the above plot I did not understand that the values in data_kf_0001.txt
were being interpolated in a linear ramp  so the background slope in the plot above is due to the drum slowing down, which I didnāt intend.
Hereās another run, with the data_kf_0001.txt
set up as follows:
5, 1
linear
! time var1
0.0 0
0.5 0
1.0 20.0
6.0 20.0
6.5 0
This keeps the drum from rotating for the first 0.5 s, while the particles settle. Having the drum stationary while the particles fall reduces the size of the initial oscillation. The AoR is just below 0.
Then from 0.5s to 1s the drum accelerates, then keeps steady for 5s. In the plot the region from 1 to 5s shows a fairly stable dynamic angle of repose, after the initial oscillations die out. The value is around 35.5 degrees as can be seen in the āDetailā plot.
From 6s to 6.5s in the plot the drum is decelerating, then stops, where the AoR reaches its final steady value of 19.5.
Another change from the previous run is that I fixed an error in the code, it is possible that the des_pos_new
array is larger than the number of particles, so my computation of the mean position was subject to extra junk at the end of the array. I changed my usr1_des.f
to this:
mean = sum(des_pos_new(1:particles,:), 1) / particles
write(*,*) "ANGLE", 180*atan2(mean(2), mean(1))/PI + 90
Very nice Charles!
The second line of the keyframe file sets the interpolation type between data points. Available options are linear
for linear interpolation and step
for no interpolation. The Conveyor_belts.pdf
file in the conveyor dem tutorial has more details.
This animation shows the computed (statistical) angle of repose as a series of parallel lines drawn across the image. Note that an angle is always computed, even when the particles are not in a settled state, hence the strange result in the first few frames.
When the particles are settled, the fit looks quite good.
When the top surface is not flat, itās a little harder to decide.
It will be up to you to determine whether this approach is adequate for your purposes.
Charles,
How did you create the frames for this animation? The lines superimposed on the drum are a really nice visual aid!
Thanks Julia.
It was a bit adhoc, I saved the output as an āimage stackā, then used a little Python script to draw the lines using the ImageMagick utility. Note that there is a library for calling ImageMagick from Python but I didnāt use it, my script just computed the line endpoint coordinates and issued commands like:
convert draw "line 0,0 1086,858" draw "line 1086,0 0,858" draw "line 543,0 543,858" x.png y.png
using os.system
calls (primitive but effective). Then, after creating all the annotated frames, I put them together using ffmpeg
, with a command like:
ffmpeg r 30 f image2 i foo%03d.png vcodec libx264 pix_fmt yuv420p z.mp4
Iām sure there are many better ways to do this, but this was my adhoc approach using the utilities that were already installed and familiar to me.
ā Charles
Charles,
I see, so I assume some component of your Python script takes the drumās geometry as an input (diameter/circular intersection with your lines) to calculate the line endpoint coordinates?
I havenāt heard of the ImageMagick utility before so thank you for bringing it up!
Exactly. In this case it was pretty simple because the center of the drum is right in the center of the image, and I measured the radius empirically (in pixels). Iād share the script I used but Iām afraid I didnāt save it after I was done!
Charles,
No worries, I was mostly interested in your thought process/workflow overall, which makes sense now!