spheres pass facet

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

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

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.

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

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

> 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

Thank you Jan!.