← Back to team overview

yade-users team mailing list archive

[Question #293955]: How to generate dense packing of aggregates

 

New question #293955 on Yade:
https://answers.launchpad.net/yade/+question/293955

Hello everyone,
I am generating a pack of spheres (lets call them aggregates), using makeCloud. Each of those spheres consists of a packing of smaller particles glued together using cohesive material model. I am using randomDensePack to generate this packing of particles within each sphere.

Now, how can I get an overall dense and mechanically stable packing of this system ? Going through some previous questions, I understand that I might have to use some type of compaction or triaxial controller, but I am not sure how it will fit in the overall flow of my simulation. Here is my minimal working example script:

 
from yade import pack

#############################################################################
# Set up run
#############################################################################
run_name="PeriodicCohesive_MWE"
data_root_dir="."


#############################################################################
# Materials
#############################################################################
plate_material=FrictMat(
    young=200e9
   ,poisson=0.3
   ,density=8000
   ,frictionAngle=radians(30)
   ,label='plate_mat')
O.materials.append(plate_material)


sample_material=CohFrictMat(
    young=4e9
   ,poisson=0.25
   ,density=1400
   ,frictionAngle=radians(30)
   ,normalCohesion=1e8*1.2*1e10
   ,shearCohesion=.4e8*1.2*1e10
   ,momentRotationLaw=True
   ,label='sample_mat')
O.materials.append(sample_material)



#############################################################################
# Component dimensions and operating condition
#############################################################################
# Coal dimension
sample_diameter=0.825e-3
sample_radius=sample_diameter/2.0
# Sub-particle dimension
particle_diameter=200.e-6
particle_radius=particle_diameter/2.

############################################################################# 
# grinding plate dimension
#############################################################################

rotvel=2./3.*pi*(1.5+0.5)*.254


#############################################################################
# Periodic Geometry
#############################################################################


# Set up periodic boundary conditions
O.periodic=True
xExt=7*sample_diameter
yExt=5.*sample_diameter*2 #to block the periodicity in y direction
zExt=xExt
xLim=xExt
yLim=yExt/4
zLim=zExt
O.cell.hSize=Matrix3(
  xExt, 0, 0,
  0, yExt, 0,
  0, 0, zExt)


length=xExt
height=yExt
width=zExt

# Top and bottom plate thickness
thickness=0.1*height



#########
# Bottom
#########

bottomBoxes = []
for ix in (0,1,2):
 for iz in (0,1,2):
   bottomBoxes.append(box( # create 3x3 boxes with 1/3 cell size
   center=(xExt/6.*(1+2*ix),yLim - thickness/2.0,zExt/6.*(1+2*iz))
   ,extents=(xExt/6.,thickness/2.0,zExt/6.)
   ,wire=False
   ,material='plate_mat'
   ))

bottom_id,bottom_ids = O.bodies.appendClumped(bottomBoxes) # bottom_id is the clump id,


O.bodies[bottom_id].state.blockedDOFs='xyzXYZ'


#############################################################################
# Particle Packing
#############################################################################

min_corner= (0,yLim,0)
max_corner= (xLim, yExt-yLim, zLim)

sp=pack.SpherePack()
sp.makeCloud( min_corner,max_corner, psdSizes = [500e-6,600e-6,700e-6,800e-6,1000e-6,1200e-6,1400e-6,], 
psdCumm=[0.0977,0.1,0.12,0.18,0.39,0.68,1.0], periodic=False, seed =1)

print "Generated ",len(sp)," aggregates"

###########################################
# Sample
###########################################
for s in sp:
  sphere=pack.inSphere((s[0][0],s[0][1],s[0][2]),s[1])
  sp1=pack.randomDensePack(
      sphere
     ,spheresInCell=5000
     ,radius=particle_radius
     ,memoizeDb='/tmp/bigParticlesCache.sqlite'
     ,returnSpherePack=True
     )

  sp1.toSimulation(material='sample_mat',color=(0.9,0.8,0.6))
  print 'Generated ',len(sp1),' particles'

Gl1_Sphere(stripes=True)


#########
 # Top
#########
topBoxes = []
for ix in (0,1,2):
 for iz in (0,1,2):
   topBoxes.append(box( # create 3x3 boxes with 1/3 cell size
   center=(xExt/6.*(1+2*ix),yExt - yLim + thickness/2.0,zExt/6.*(1+2*iz))
   ,extents=(xExt/6.,thickness/2.0,zExt/6.)
   ,wire=False
   ,material='plate_mat'
   ))

top_id,top_ids = O.bodies.appendClumped(topBoxes) # top_id is the clump id,


O.bodies[top_id].state.blockedDOFs='xzXYZ'



one_ball_downforce= -15.0 # kg total

O.forces.addF(top_id,(0,one_ball_downforce,0),permanent=True)


#############################################################################
# Run the simulation
#############################################################################
r_int = 1.1 # interaction radius

itercheck=10000
iterwrite=10000
iterforceCheck = 100000
iterStart=0 
tStart=0.


O.dt=0.5*PWaveTimeStep()

O.engines=[
  ForceResetter(),
  InsertionSortCollider([
    Bo1_Sphere_Aabb(aabbEnlargeFactor=r_int, label='bo1s')
   ,Bo1_Box_Aabb()
   ], allowBiggerThanPeriod=True
   ),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=r_int, label='ig2ss')
    ,Ig2_Box_Sphere_ScGeom6D() 
    ],
    [Ip2_FrictMat_FrictMat_FrictPhys()
     ,Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,label="cohesiveIp")
    ],
    [Law2_ScGeom_FrictPhys_CundallStrack()
     ,Law2_ScGeom6D_CohFrictPhys_CohesionMoment()
    ]
  ), # End InteractionLoop 
  NewtonIntegrator(damping=0.8,gravity=(0.,-9.81,0.)),

  PyRunner(command='forceCheck()', iterPeriod=iterforceCheck, label='forceCheckEngine'),
  

]

O.step()



def forceCheck():
  
  reactionForce = sum( ( O.forces.f(i) for i in top_ids ), Vector3.Zero )
  print reactionForce
  O.bodies[top_id].state.vel[0]= rotvel
  forceCheckEngine.dead = True




for i in O.interactions: 
    if not isinstance(i.phys,CohFrictPhys): # ... only ones with CohFrictPhys interactions
			continue
    i.phys.normalAdhesion = 1.2
    i.phys.shearAdhesion  = 0.48
    i.phys.unp = 0 #to make it initially force free	
    

O.run()




-- 
You received this question notification because your team yade-users is
an answer contact for Yade.