New question #695929 on Yade:
https://answers.launchpad.net/yade/+question/695929
Hi all,
I would like to simulate the compaction of granular material. I am using
yadedaily and my code is shown below. I have 2 questions:
1- I see that smaller spheres escape from the compaction mold during the
compaction. How to prevent any sphere from moving out of the mold (i.e. facets)?
2- I want to simulate packing at high compaction force. In my current code,
after around 6000 iteration, the packing "explode" and spheres get out of the
modl. Why this happens? Is it a numerical issue or this is what suppose to
happen based on the Cundall-Strack contact model?
I appreciate your help and comments.
Thank you,
Othman
--------------------------------------------------------
from yade import pack, export, ymport, plot
from pprint import pprint
import numpy as np
import matplotlib.pyplot as plt
import time
tic=time.time()
O.materials.append(FrictMat(young = 1e6, poisson = 0.45,frictionAngle = 0.349,
density=2340))
SG=2.34
##cylinder dimensions
radiuscyl=(500e-6/2)
heightcyl=610e-6
##center of cylinder
cx=0
cy=0
cz=0
##Initial cube dimensions ###
mnx=cx-(radiuscyl*1.1)
mny=cy-(radiuscyl*1.1)
mnz=0
mxx=cx+(radiuscyl*1.1)
mxy=cy+(radiuscyl*1.1)
mxz=heightcyl*1.1
############################ spheres #############################
sp=pack.SpherePack()
##### sizes and distribution are from gradation curve of basalt aggregates
#######
sizes=1e-6*np.array([1.58*0.98,1.58,3.63,9.16,15.76,22.94,25]) #Diameters of
portlandite
passing=[0,0.1,0.25,0.5,0.75,0.9,1]
sp.makeCloud((mnx,mny,mnz),(mxx,mxy,mxz),psdSizes=sizes,psdCumm=passing)
#### cylinder extraction
pred=pack.inCylinder((cx,cy,cz),(cx,cy,heightcyl),radiuscyl)
spFilter=filterSpherePack(pred,sp,Material=Material, returnSpherePack=True)
print (len (spFilter))
spFilter.toSimulation(color=(0.533, 0.803, 0.780))
print ("runtime = ", time.time()-tic)
mass=utils.getSpheresMass()
############################ facets #############################
facets=geom.facetCylinder((cx,cy,heightcyl/2),radiuscyl,heightcyl,segmentsNumber=150,wallMask=4)
cylinder=O.bodies.append(facets)
yade.qt.View()
##creating disks
d1=geom.facetCylinder((cx,cy,heightcyl),radiuscyl*0.99,0,segmentsNumber=300,wallMask=1)
d2=geom.facetCylinder((cy,cx,cz),radiuscyl*0.99,0,segmentsNumber=300,wallMask=1)
disk1IDs= O.bodies.append(d1)
disk2IDs= O.bodies.append(d2)
for i in disk1IDs:
body= O.bodies[i]
body.state.vel = (0,0,-3.5)
for n in disk2IDs:
body= O.bodies[n]
body.state.vel = (0,0,0)
############################# compaction #############################
O.dt=.5*utils.PWaveTimeStep()
O.engines=[
ForceResetter(),
InsertionSortCollider([
Bo1_Sphere_Aabb(),
Bo1_Facet_Aabb()
]),
InteractionLoop(
[
Ig2_Sphere_Sphere_ScGeom(),
Ig2_Facet_Sphere_ScGeom(),
],
[
Ip2_FrictMat_FrictMat_FrictPhys(),
Ip2_FrictMat_FrictMat_FrictPhys(),
],
[
Law2_ScGeom_FrictPhys_CundallStrack(),
],
),
NewtonIntegrator(damping=.3),
# VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=500),
PyRunner(iterPeriod=500,command='force()',initRun=True),
]
O.run()
stress=[]
Thickness=[]
Packing_d=[]
def force():
f1= [O.forces.f(i)[2] for i in disk1IDs]
f=np.mean(f1)
s=f/(pi*radiuscyl**2) #stress N/m2
stress.append(s)
thickness=(O.bodies[disk1IDs[1]].state.pos[2])-(O.bodies[disk2IDs[1]].state.pos[2])
packing_density=mass/(thickness*pi*radiuscyl**2)/997/SG
print("stress, thickness, packing density ",s,thickness,packing_density)
Thickness.append(thickness)
Packing_d.append(packing_density)
plot.addData(applied_stress=s,thickness=thickness,packing_density=packing_density)
plot.plots={('applied_stress'):('packing_density')}
plot.plot()
#np.savetxt('compaction_data.txt',np.transpose([stress,Packing_d,Thickness]),delimiter=',')
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
