yade-users team mailing list archive

[Othman Sh <question695929@xxxxxxxxxxxxxxxxxxxxx>]

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.