Data Files for “Rolling Resistance Linear Contact Model: Repose Angle” Problem
PFC2D Data Files
make_system.p2dat
;===========================================================================
; fname: make_system.p2dat
;===========================================================================
model new
model title 'Rolling Resistance Linear Contact Model: Repose angle'
model random 10001
; container and balls attributes
[cyl_rad = 50.0e-3 ] ; container half-width [m]
[cyl_height = 100.0e-3] ; container height [m]
[ravg = 1.5e-3 ] ; ball average radius [m]
[rdev = 0.5 ] ; ball relative radius deviation [-]
[dens = 2.0e3 ] ; ball density [kg/m3]
[poros = 0.3 ] ; initial porosity [-]
; contact model properties
[emod = 1.0e6] ; Young's modulus (deformability method) [Pa]
[kratio = 2.0 ] ; stiffness ratio (deformability method) [-]
[fric = 0.01 ] ; friction coefficient [-]
[dpnr = 0.2 ] ; normal critical damping ratio [-]
[dpsr = 0.2 ] ; shear critical damping ratio [-]
[dpm = 3 ] ; dashpot mode [-]
[rfric = 0.0 ] ; rolling resistance coefficient [-]
; set domain extent and boundary conditions
model domain extent [-4.0*cyl_rad] [4.0*cyl_rad] ...
[-4.0*ravg] [2.0*cyl_height] ...
condition destroy
; set the CMAT default slots
contact cmat default model rrlinear method deformability ...
emod @emod ...
kratio @kratio ...
property fric @fric ...
dp_nratio @dpnr ...
dp_sratio @dpsr ...
dp_mode @dpm ...
rr_fric @rfric
; generate the base plane and the container
wall generate id 1 ...
name 'base_plane' ...
plane
wall create id 100 ...
name 'container' ...
vertices [-cyl_rad] 0.0 [-cyl_rad] @cyl_height ...
[-cyl_rad] @cyl_height [ cyl_rad] @cyl_height ...
[cyl_rad] @cyl_height [ cyl_rad] 0.0
; distribute balls inside the container region
ball distribute porosity @poros ...
resolution @ravg ...
radius [1.0-0.5*rdev] [1.0+0.5*rdev] ...
box [-1.0*cyl_rad] [1.0*cyl_rad] ...
0.0 @cyl_height
; bring to rest under gravity loading
ball attribute density @dens damp 0.7
model cycle 1000 calm 10
model gravity 9.81
model solve
ball attribute damp 0.0
model save 'system_ini.p2sav'
program return
;===========================================================================
; eof: make_system.p2dat
move_container.p2dat
;===========================================================================
; fname: move_container.p2dat
;===========================================================================
[cyl_vel = 5.0e-2] ; upwards velocity [m/s]
model mechanical time-total 0.0
contact cmat default property fric @fric rr_fric @rfric
contact property fric @fric rr_fric @rfric
wall attribute velocity-y @cyl_vel range id 100 by wall
ball attribute displacement multiply 0.0
model solve time [0.75*cyl_height/cyl_vel]
wall attribute velocity-y 0.0 range id 100 by wall
model solve ratio-average 1e-5
model save [string.build('system_final-fric%1-rfric%2.p2sav',fric,rfric)]
program return
;===========================================================================
; eof: move_container.p2dat
do_analysis.p2dat
;===========================================================================
; fname: do_analysis.p2dat
;===========================================================================
model restore 'system_ini'
[fric = 0.5]
[rfric = 0.1]
program call 'move_container.p2dat'
;
model restore 'system_ini'
[fric = 0.5]
[rfric = 0.5]
program call 'move_container.p2dat'
;
model restore 'system_ini'
[fric = 0.5]
[rfric = 1.0]
program call 'move_container.p2dat'
;===========================================================================
; eof: do_analysis.p2dat
do_analysis.py
#===========================================================================
# do_analysis.py
#===========================================================================
import itasca
fric = [0.25,0.5]
rfric = [0.05,0.1,0.25,0.5,0.6]
for f in fric:
for rf in rfric:
itasca.command("""
model restore 'system_ini'
[fric = {0}]
[rfric = {1}]
program call 'move_container.p2dat'
""".format(f,rf)
)
#===========================================================================
# eof: do_analysis.py
PFC3D Data Files
make_system.p3dat
;=============================================================================
; fname: make_system.p3dat
;=============================================================================
model new
model title 'Rolling Resistance Linear Contact Model: Repose angle'
model random 10001
; container and balls attributes
[cyl_rad = 50.0e-3 ] ; container radius [m]
[cyl_height = 100.0e-3] ; container height [m]
[ravg = 3.0e-3 ] ; ball average radius [m]
[rdev = 0.25 ] ; ball relative radius deviation [-]
[dens = 2.0e3 ] ; ball density [kg/m3]
[poros = 0.45 ] ; initial porosity [-]
; contact model properties
[emod = 1.0e6]
[kratio = 2.0 ]
[fric = 0.05 ]
[dpnr = 0.2 ]
[dpsr = 0.2 ]
[dpm = 3 ]
[rfric = 0.0 ]
; set domain extent and boundary conditions
model domain extent [-4.0*cyl_rad] [4.0*cyl_rad] ...
[-4.0*cyl_rad] [4.0*cyl_rad] ...
[-4.0*ravg] [2.0*cyl_height] ...
condition destroy
; set the CMAT default slots
contact cmat default model rrlinear method deformability ...
emod @emod ...
kratio @kratio ...
property fric @fric ...
dp_nratio @dpnr ...
dp_sratio @dpsr ...
dp_mode @dpm ...
rr_fric @rfric
; generate the base plane and the container
wall generate id 1 ...
name 'base_plane' ...
plane
wall generate id 100 ...
name 'container' ...
cylinder ...
base (0.0,0.0,0.0) ...
axis (0.0,0.0,0.1) ...
height @cyl_height ...
radius @cyl_rad ...
cap false true ...
one-wall
; distribute balls inside the container region
ball distribute porosity @poros ...
resolution @ravg ...
radius [1.0-0.5*rdev] [1.0+0.5*rdev] ...
box [-1.0*cyl_rad] [1.0*cyl_rad] ...
[-1.0*cyl_rad] [1.0*cyl_rad] ...
0.0 @cyl_height ...
range cylinder end-1 (0.0,0.0,0.0) ...
end-2 (0.0,0.0,@cyl_height) ...
radius @cyl_rad ...
extent
; bring to rest under gravity loading
ball attribute density @dens damp 0.7
model cycle 1000 calm 10
model gravity 9.81
model solve
ball attribute damp 0.0
model save 'system_ini'
return
;=============================================================================
; eof: make_system.p3dat
move_container.p2dat
;=============================================================================
; fname: move_container.p3dat
;=============================================================================
[cyl_vel = 5.0e-2]
model mechanical time-total 0.0
contact cmat default property fric @fric rr_fric @rfric
contact property fric @fric rr_fric @rfric
wall attribute velocity-z @cyl_vel range id 100 by wall
ball attribute displacement multiply 0.0
model solve time [0.6*cyl_height/cyl_vel]
wall attribute velocity-z 0.0 range id 100 by wall
model solve ratio-average 1e-4
model save [string.build('system_final-fric%1-rfric%2.p3sav',fric,rfric)]
return
;=============================================================================
; eof: move_container.p3dat
do_analysis.py
import itasca
fric = [0.2,0.6]
rfric = [0.1,0.2,0.4,0.6,0.8]
for f in fric:
for rf in rfric:
itasca.command("""
model restore \'system_ini\'
[fric = {0}]
[rfric = {1}]
program call \'move_container\' suppress
""".format(f,rf)
)
post_treat.py
import utils
import matplotlib.pyplot as plt
fric = [0.2,0.6]
rfric = [0.1,0.2,0.4,0.6,0.8]
val = []
conf = []
leg = []
for f in fric:
ra = []
ci = []
for rf in rfric:
p,c = utils.fit_cone('system_final-fric{0}-rfric{1}.p3sav'.format(f,rf))
ra += [p]
ci += [c]
val += [ra]
conf += [ci]
leg += [r'$\mu={0}$'.format(f)]
plt.close('all')
for i, v in enumerate(val):
plt.errorbar(rfric , val[i],yerr=ci[i])
plt.legend(leg,loc=2)
plt.title("Repose Angle vs Rolling Friction Coefficient")
plt.xlabel("Rolling Friction Coef. [-]")
plt.ylabel("Angle of Repose [$^{\circ}$]")
plt.axis([0, 0.9, 0, 40])
plt.grid(True)
plt.draw()
utils.py
import itasca
from itasca import ballarray as ba
import numpy as np
from scipy.spatial.kdtree import KDTree
from scipy import optimize
from scipy import stats
import os
def cone_eq(data,a,b):
return a*(1-np.sqrt(data[:,0]**2+data[:,1]**2)/
(a*np.tan((90.0-b)*np.pi/180.0)))
def locally_extreme_points(coords, data, neighbourhood,
lookfor = 'max', p_norm = 2.):
'''
Find local maxima of points in a pointcloud.
Ties result in both points passing through the filter.
Not to be used for high-dimensional data. It will be slow.
coords: A shape (n_points, n_dims) array of point locations
data: A shape (n_points, ) vector of point values
neighbourhood: The (scalar) size of the neighbourhood in which to search.
lookfor: Either 'max', or 'min', depending on whether you
want local maxima or minima
p_norm: The p-norm to use for measuring distance
(e.g. 1=Manhattan, 2=Euclidian)
returns
filtered_coords: The coordinates of locally extreme points
filtered_data: The values of these points
'''
assert coords.shape[0] == data.shape[0], \
'You must have one coordinate per data point'
extreme_fcn = {'min': np.min, 'max': np.max}[lookfor]
kdtree = KDTree(coords)
neighbours = kdtree.query_ball_tree(kdtree,
r=neighbourhood, p = p_norm)
i_am_extreme = [data[i]==extreme_fcn(data[n])
for i, n in enumerate(neighbours)]
extrema, = np.nonzero(i_am_extreme) # This line just saves time on indexing
return extrema,coords[extrema], data[extrema]
def fit_cone(fname):
'''
Fit a cone
fname: A save file name
returns
t_fit: The repose angle
ci: The student-t 95% confidence interval
'''
itasca.command("""
model restore \'{}\'
""".format(fname)
)
bname = os.path.splitext(fname)[0]
sname = bname + '-ps.p3sav'
positions = ba.pos()
radii = ba.radius()
rmin= np.amin(ba.radius())
rmax= np.amax(ba.radius())
ravg= 0.5*(rmin+rmax)
x,y,z = positions.T
zmax = z + radii
coords = np.vstack((x, y)).T
extrema, newcoords,val = locally_extreme_points(coords, z,
2.0*ravg,
lookfor = 'max',
p_norm = 2.)
damp = ba.damp()
damp[extrema] = 1
ba.set_damp(damp)
xc = x[extrema]
yc = y[extrema]
zc = zmax[extrema]
# Fit a cone on the local extrema.
# Use scipy.optimize.curve_fit to estimate the confidence interval
cC = np.amax(zc)
mask_up = (cC - 2.0*rmax) > zc
mask_lo = zc > (2.0*rmax)
xm = xc[mask_up & mask_lo]
ym = yc[mask_up & mask_lo]
zm = zc[mask_up & mask_lo]
data = np.vstack((xm, ym)).T
p0 = [cC,40.0] # initial guess for optimization
params, pcov = optimize.curve_fit(cone_eq, data, zm, p0)
# student-t value for the dof and confidence level
n = len(zm) # number of data points
p = len(params) # number of parameters
dof = max(0, n - p) # number of degrees of freedom
alpha = 0.05 # 95% confidence interval = 100*(1-alpha)
tval = stats.distributions.t.ppf(1.0-alpha/2., dof)
ci = []
for i, p,var in zip(range(n), params, np.diag(pcov)):
sigma = var**0.5
ci.append(sigma*tval)
print('p{0}: {1} [{2} {3}]'.format(i, p,
p - sigma*tval,
p + sigma*tval))
c_fit, t_fit = params
R_fit = c_fit*np.tan((90.0-t_fit)*np.pi/180.0)
itasca.command(
"""
ball group \'mask\' range position-z {zmin} {zmax} not
[repose_angle = {ra}]
geometry generate cone base (0,0,0) ...
axis (0,0,1) ...
height {c} ...
radius {r} ...
cap false
model save \'{s}\'
""".format(ra=t_fit,c=c_fit,r=R_fit,
zmin = 2.0*rmax,
zmax = (cC - 2.0*rmax),s=sname)
)
return [t_fit,ci[1]]
⇐ Rolling Resistance Linear Contact Model: Repose Angle | Adhesive Rolling Resistance Linear Contact Model: Repose Angle ⇒
| Was this helpful? ... | PFC 6.0 © 2019, Itasca | Updated: Nov 19, 2021 |