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
cylinder = geom.facetCylinder(center=center, radius=diameter/2, height=height, segmentsNumber=80, wallMask=6)

# add cylinder to simulation
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)

# add funnel to simulation
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
cylinder = geom.facetCylinder(center=center, radius=diameter/2, height=height, segmentsNumber=80, wallMask=6)

# add cylinder to simulation
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)

# add funnel to simulation
O.bodies.append(funnel)

# add sphere packing
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')

# add properties
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
   if body.shape.radius == 0.01575/2 :
       body.shape.color = (0,0,1) #blue
       body.material = gravel
   if body.shape.radius == 0.011/2:
       body.shape.color = (1,0,0) #red
       body.material = gravel
   if body.shape.radius == 0.007125/2:
       body.shape.color = (0,1,0) #green
       body.material = gravel
   if body.shape.radius == 0.003555/2:
       body.shape.color = (1,1,0) #yellow
       body.material = gravel
   if body.shape.radius == 0.00160/2 :
       body.shape.color = (1,0,1) #magenta
       body.material = gravel
   if body.shape.radius == 0.0008/2 :
       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
        PyRunner(command='addPlotData()', iterPeriod=100)
]
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
def addPlotData():
 # 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
 plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)

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