yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #15330
[Question #649832]: Porosity and friction angle
New question #649832 on Yade:
https://answers.launchpad.net/yade/+question/649832
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()
--
You received this question notification because your team yade-users is
an answer contact for Yade.