# I try to pack the sphere in the cylinder by gravity deposition. But it stop when it touches the funnel.

Hello,
I create a sphere packing position in the 1st script. However in second script, I import the position from the 1st script. I reduce the sphere value. I am confuse on why once the sphere hit the funnel the simulation stop right away.

### 1st script ###
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 = 0.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.60 mm 51x1400 = 71,400 particles
n5 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.0016/2, num=1400)

# 0.8 mm 51x4140 = 211,140 particles
n6 = sp.makeCloud(minCorner=corner[0], maxCorner=corner[1], rMean=0.0008/2, num=4140)

sp.toSimulation()

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

### 2nd script ###
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 = 27000, label = 'gravel')
asphalt_binder = CohFrictMat(young = 1e7, poisson = 0.25, density = 10600, 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 = (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, -981), 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 = 0.1 * 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() < 0.01:
O.pause()
export.text("gra_depo_pos.txt")
plot.saveDataTxt('bbb.txt')
# 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:
For:
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Huan (huan-liu) said on 2023-05-26: #1

I did modify the 1st script for lesser sphere to see if the same reaction happens. Around 60k sphere particles also have the same reaction once the sphere hit the funnel the simulation automatically.

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

Most likely option is that it met the condition " unbalancedForce() < 0.01" in checkUnbalanced() function, where the simulation is paused.
You can print the value before O.pause()
Cheers
Jan

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

>What value to print?
Do you mean print unbalancedForce before O.pause(). And what would printing those value does.

> But how do I prevent the sphere from stopping when it touches the funnel or how do I let it complete the simulation for me to get the position when the sphere is static within the cylinder?

 Revision history for this message Jan Stránský (honzik) said 10 hours ago: #4

>> You can print the value before O.pause()
>What value to print?
> Do you mean print unbalancedForce before O.pause(). And what would printing those value does.

Printing that value will:
- give you evidence / proof that the problem is really the unbalanced force condition (print just before O.pause)
- just print the unbalanced force value, so you can compare it to the condition limit and adjust it accordingly

> But how do I prevent the sphere from stopping when it touches the funnel or how do I let it complete the simulation for me to get the position when the sphere is static within the cylinder?

Obviously (if the unbalancedForce is really the problem), choose a different value or different condition.
Or add another condition, e.g. that most of the particles are already "low enough".
Or ... really many options here depending on actual needs, requirements, resources etc.

Cheers
Jan