Sphere packing location

Hello,

I create two python script.
1. Sphere Packing Position Saving
2. Importing the sphere packing from 1. to do gravity deposition
>>The Problem I encounter is that I try to move the ( zi = minz + l * gridsize ) lower by dividing it by 8. However, after I import the postion in 2nd python script the sphere packing burst out which I am confuse on what is making them burst apart.

### SCRIPT 1 ###
import random
import math
from yade import geom, pack, utils, plot, ymport, export
import numpy as np

# Define cylinder with funnel parameters
center = (0, 0, 0)
diameter = 0.102
height = 0.18

# create cylindrical body with radius 0.102 m and height 0.064 m

O.bodies.append(cylinder)

# Define cylinder with funnel parameters
center1 = (0,0,height/2)
dBunker = 0.4
dOutput = 0.102
hBunker = 0
hOutput = 0.15
hPipe = 0

# create funnel as a bunker with diameter 0.102 m, height 0.064 m
funnel = geom.facetBunker(center=center1, dBunker=dBunker, dOutput=dOutput, hBunker=hBunker,hOutput=hOutput, hPipe=hPipe, segmentsNumber=80, wallMask=4)

O.bodies.append(funnel)

def makeRandomBlocks(diameter, gridsize, minz, numblocks):
bunkerLimit = diameter/2
#
# Make grid blocks inside a bunker
#
grid = np.arange(-bunkerLimit, bunkerLimit, gridsize)
inGrid = []
for yi in grid:
for xi in grid:
c = [[xi, yi], [xi+gridsize, yi], [xi+gridsize, yi+gridsize], [xi, yi+gridsize]]
#
# Check if the block corners are outside of the bunker or not
#
out = False
for p in c:
r = (p[0]**2 + p[1]**2)**0.5
if r > bunkerLimit:
out = True
#
# If the block is inside the bunker, keep it
#
if not out:
inGrid.append(c)
layer = math.ceil(numblocks / len(inGrid))
blocks = []
for l in range(layer):
for g in inGrid:
zi = minz + l*gridsize
p1 = g[0].copy()
p1.append(zi)
p2 = g[2].copy()
p2.append(zi+gridsize)
blocks.append([p1, p2])
random.shuffle(blocks)
return blocks

minZ = 2.35
dbunker = 0.4
gridSize = dbunker/8
numblocks = 51
blockList = makeRandomBlocks(dbunker, gridSize, minZ, numblocks)

for i in range(numblocks):
print("Making cloud block", i+1, "/", numblocks)
corner = blockList.pop()
sp = pack.SpherePack()

# 15.75 mm 13 particles
if (i < 13):
n1 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.01575/2, num=1)

# 11 mm 51 particles
n2 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.011/2, num=1)

# 7.125 mm 51x11 = 561 particles
n3 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.007125/2, num=11)

# 3.555 mm 51x100 = 5,100 particles
n4 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.003555/2, num=100)

# 1.77 mm 51x360 = 18,360 particles
n5 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.00177/2, num=360)

# 0.6 mm 51x19000 = 969,000 particles
n6 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.0006/2, num=19000)

# 0.3 mm 51x78510 = 4,004,010 particles
n7 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.0003/2, num=78510)

sp.toSimulation()

export.text("testCloud.txt")
########################################################

### Script 2 ###
import random
import math
from yade import geom, pack, utils, plot, ymport, export
import numpy as np

# Define cylinder with funnel parameters
center = (0, 0, 0)
diameter = 0.102
height = 0.18

# create cylindrical body with radius 0.102 m and height 0.064 m

O.bodies.append(cylinder)

# Define cylinder with funnel parameters
center1 = (0,0,height/2)
dBunker = 0.4
dOutput = 0.102
hBunker = 0
hOutput = 0.15
hPipe = 0

# create funnel as a bunker with diameter 0.102 m, height 0.064 m
funnel = geom.facetBunker(center=center1, dBunker=dBunker, dOutput=dOutput, hBunker=hBunker,hOutput=hOutput, hPipe=hPipe, segmentsNumber=80, wallMask=4)

O.bodies.append(funnel)

O.bodies.append(ymport.textExt('testCloud.txt',format='x_y_z_r'))

# materials Properties
gravel = CohFrictMat(young = 1e7, poisson = 0.25, density = 2700, label = 'gravel')
asphalt_binder = CohFrictMat(young = 1e7, poisson = 0.25, density = 1060, frictionAngle = radians(40), normalCohesion = 2e5, shearCohesion = 2e5, label = 'asphalt_binder')

O.materials.append(gravel)
O.materials.append(asphalt_binder)

# give color and properties to shpere
for body in O.bodies:
if not isinstance(body.shape, Sphere):
continue
body.shape.color = (0,0,1) #blue
body.material = gravel
body.shape.color = (1,0,0) #red
body.material = gravel
body.shape.color = (0,1,0) #green
body.material = gravel
body.shape.color = (1,1,0) #yellow
body.material = gravel
body.shape.color = (1,0,1) #magenta
body.material = gravel
body.shape.color = (1,1,1) #white
body.material = gravel
body.shape.color = (0,0,0) #black
body.material = asphalt_binder

O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
InteractionLoop(
# handle sphere+sphere and facet+sphere collisions
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.4),
# call the checkUnbalanced function (defined below) every 2 seconds
PyRunner(command='checkUnbalanced()', realPeriod=2),
# call the addPlotData function every 200 steps
]
O.dt = .5 * PWaveTimeStep()

# enable energy tracking; any simulation parts supporting it
# can create and update arbitrary energy types, which can be
# accessed as O.energy['energyName'] subsequently
O.trackEnergy = True

# if the unbalanced forces goes below .05, the packing
# is considered stabilized, therefore we stop collected
# data history and stop
def checkUnbalanced():
if unbalancedForce() < .05:
O.pause()
plot.saveDataTxt('bbb.txt.bz2')
# plot.saveGnuplot('bbb') is also possible

# collect history of data which will be plotted
# each item is given a names, by which it can be the unsed in plot.plots
# the **O.energy converts dictionary-like O.energy to plot.addData arguments

# define how to plot data: 'i' (step number) on the x-axis, unbalanced force
# on the left y-axis, all energies on the right y-axis
# (O.energy.keys is function which will be called to get all defined energies)
# None separates left and right y-axis
plot.plots = {'i': ('unbalanced', None, O.energy.keys)}

# show the plot on the screen, and update while the simulation runs
plot.plot()

O.saveTmp()
##############################################

Question information

Language:
English Edit question
Status:
Open
For:
Assignee:
No assignee Edit question
Last query:

This question was reopened

• by Huan
 Revision history for this message Jan Stránský (honzik) said on 2023-05-25: #1

Hello,

please provide a MWE [1], M = minimal
Is 5 millions spheres needed to demonstrate your problem?

"bursting" might be caused by particle overlap.
Which coordinate "dividing it by 8" supports (if I understand it correctly)..

Cheers
Jan

 Revision history for this message Huan (huan-liu) said on 2023-05-25 #2

Thanks Jan,
I was able to modify it by lowering the min z.
I think 5 millions spheres is not relevant, but my prof. want all of it even though my pc aren't able to generate

 Revision history for this message Huan (huan-liu) said on 2023-05-25: #3

Hey,

Jan out of curiosity if the sphere that are generated are so small for example my asphalt binder is 0.0003 wouldn't it just go through the funnel and the cylinder for gravity deposition part?

 Revision history for this message Jan Stránský (honzik) said on 2023-05-25: #4

I understand that for actual simulation, you use that many particles.
My point was that for this forum and for higher probability of good answer (I was not able to run the script on my laptop, so I was only guessing), always try to create a script focusing just on the problem, trying to make it minimal (both from the point of view of code and resources needed for run).
As a byproduct, by creating such MWE, often you can find the problem yourself.

"Small" particles (in fact any particle) can "go through" a barrier, it depends on many factors, mainly:
- if velocity of the particle is not as high as that in one time step it fully goes from one side of the barrier to the other, without any possibility of contact detection
- if stiffness / force / inertia interplay is not such that the barrier cannot prevent the particle from passing through.

Cheers
Jan

 Revision history for this message Gem Herrera (gem24) said on 2023-05-26: #5

Suppose velocity is not as high as that in one time step, it fully goes from one side of the barrier to the other, without of contact detection
And inertia interplay is not such that the barrier cannot stop the particle from passing through.
Regards
https://proassemblage.com/