spheres pass facet

Asked by Othman Sh on 2021-03-08

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=',')

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Othman Sh
Solved:
2021-03-09
Last query:
2021-03-09
Last reply:
2021-03-09
Jan Stránský (honzik) said : #1

Hello,
for both cases, smaller value of O.dt should help
cheers
Jan

Othman Sh (othman-sh) said : #2

Hi Jan,

I tried reducing the O.dt from O.dt=.5*utils.PWaveTimeStep() to O.dt=.005*utils.PWaveTimeStep(). Although I reduced it 100 times, I still get spheres pass through the facets. How do I know what O.dt is suitable for my case? Any other solution?

Thanks,
Othman

Robert Caulk (rcaulk) said : #3

You also would need a higher stiffness for the walls than the spheres. However, there is no way to know what stiffness you are currently using for your walls since the code you post probably has some kind of python error saying "Material" is not defined. Please send an MWE for more help.

Jan Stránský (honzik) said : #4

Another option is to use boxes instead of facets - boxes has finite thickness and can help in this situation.
Jan

Othman Sh (othman-sh) said : #5

Hi Jan and Robert,

My experimental test set-up is on cylinder compaction so I really need to use a cylinder in my simulation. I have defined "Material" and the code copied below works now.

Robert, how to change the stiffness of the walls? I am assuming I should assign a different material and elastic modulus for the facets but I checked the facetCylinder function [1] and it has no "material" parameter.

Thanks for the help,
Othman

[1] https://yade-dem.org/doc/yade.geom.html#yade.geom.facetCylinder

------------------

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()

Material=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=',')

Jan Stránský (honzik) said : #6

> how to change the stiffness of the walls? I am assuming I should assign a different material and elastic modulus for the facets

exactly. Because Yade uses harmonic mean, stiffness of walls can be safely "infinite"

> but I checked the facetCylinder function [1] and it has no "material" parameter.

material is hidden in the **kw, "(unused keyword arguments) passed to utils.facet" [1]. utils.facet [2] then have material

cheers
Jan

[2] https://yade-dem.org/doc/yade.utils.html#yade.utils.facet

Othman Sh (othman-sh) said : #7

Thank you Jan!.