# Shaking and Packing problem

I want to generate a packing from a cloud of spheres contained in a box.

Once obtained the equilibrium after the initial gravity deposition, the top of the box should go down and stop in the middle of the box. After that, I want to reverse the gravity many times in order to simulate the shaking of the spheres enclosed in the half of the box. Finally, the top facet should go down more, compress the packing, go up relaxing the spheres and, after that, I want to determine the final density/porosity of the packing.

I've already written a simple python input which i paste here. My simulation is just a draft and stops when the unbalance force is lower that a minimum value. I have basically the following questions:

1 - How can I write the condition that stops the top facet in a certain position of the box independetly from the unbalance force related to the spheres?
2 - How can I save the final configuration of the spheres which I want to load in the subsequent step with the inverted gravity?
3 - How can save/load the final position of the facet top that I want to move up in the final step of my simulation?

Here is the code:

# import yade modules that we will use below

#Materials
fr = 0.455;rho=2650.0

Mat1 = O.materials.append(ViscElMat(kn=10.0e4,ks=10.0e4,frictionAngle=fr,density=rho))
Mat2 = O.materials.append(ViscElMat(kn=10.0e4,ks=10.0e4,frictionAngle=fr,density=rho))

# create rectangular box from facets
ground=O.bodies.append(utils.wall(position=(0,0,0),axis=2,sense=0,material=Mat1))

# create empty sphere packing
sp=pack.SpherePack()

# generate spheres

sp.makeCloud((0,0,0),(10,10,20),rMean=0.5,num=1000)

# add the sphere pack to the simulation
sp.toSimulation(material=Mat2)

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
InteractionLoop(
# handle sphere+sphere and facet+sphere collisions
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
# call the checkUnbalanced function (defined below) every 2 seconds
#PyRunner(command='checkUnbalanced()',realPeriod=2),

PyRunner(command='hightcontrol()',iterPeriod=50000),

TranslationEngine(translationAxis=(0,0,-1),velocity=2,ids=top,label='trans'),

NewtonIntegrator(gravity=(0,0,-9.81),damping=0.0),

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

O.trackEnergy=True

def checkUnbalanced():
if unbalancedForce()<.05:
O.pause()

)

plot.plots={'i':('UnbalancedForce',None,'Porosity','Density',
)}

plot.plot()

O.saveTmp()

## Question information

Language:
English Edit question
Status:
Expired
For:
Assignee:
No assignee Edit question
Last query:
2020-02-06
2020-02-22
 Jan Stránský (honzik) said on 2020-02-06: #1

Hi,

welcome :-)

> 1 - How can I write the condition that stops the top facet in a certain position of the box independetly from the unbalance force related to the spheres?

similarly as checking unbalanced force:
###
def checkFacetPosition():
facet = top[0] # maybe not 0, simply one of the top facets
if facet.state.pos[2] < limit: # pos[2] is z coordinate
O.pause()
###

> 2 - How can I save the final configuration of the spheres which I want to load in the subsequent step with the inverted gravity?

export.text [1] is the easiest way
But "staying" in one script (no saving/loading) is the preferred way (although not necessary)

3 - How can save/load the final position of the facet top that I want to move up in the final step of my simulation?

also consider the one-script solution.
You can get the position of top facet as in point 1). You can save it in a file and load is by reading that file for example.

cheers
Jan

 Paolo (p4olo) said on 2020-02-06: #2

If I consider a one-script solution, how can I define a function which overwrites the gravity engine and puts +9.81?

 Jan Stránský (honzik) said on 2020-02-06: #3

use label for NewtonIntegrator:
NewtonIntegrator(...,label="integrator"),
and then you can simply do:
integrator.gravity = (0,0,+9.81)

There are several ways where and how to use that command, I cannot find any reference at the moment..

cheers
Jan

 Bruno Chareyre (bruno-chareyre) said on 2020-02-06: #4

> integrator.gravity = (0,0,+9.81)

Or, without label, O.engines[5].gravity=...
Bruno

 Paolo (p4olo) said on 2020-02-06: #5

Thank you Jan and Bruno for the hint about the NewtonIntegrator. It works. However I cannot make the top stop at a desired height. In particular I have a box which is 20 in height. I want the top to stop when it reaches 15 (so it goes down only for 5). I paste here the modified code (I made some changes from the previous code by looking at the oedometric example in the Yade Manual).
Basically I would like to reproduce the following steps:

1. Generation of the packing from a cloud.
2. Once the pack is stable (i.e. unbalancedforce<0.1) generation of the top (plate) which goes down.
3. Stop the plate at 15 height.
4. Reverse the gravity.
5. Wait the new equilibrium configuration ( i.e. unbalancedforce<0.1) and stop the simulation.

Then I need to reproduce 4 and 5 many times in order to simulate the shaking but I think that this could be easy to handle.

Here is the modified code.

# import yade modules that we will use below

#Materials
fr = 0.455;rho=2650.0
global vel
vel=2

Mat1 = O.materials.append(ViscElMat(kn=10.0e4,ks=10.0e4,frictionAngle=fr,density=rho))
Mat2 = O.materials.append(ViscElMat(kn=10.0e4,ks=10.0e4,frictionAngle=fr,density=rho))

# create rectangular box from facets
ground=O.bodies.append(utils.wall(position=(0,0,0),axis=2,sense=0,material=Mat1))

# create empty sphere packing
sp=pack.SpherePack()

# generate spheres

sp.makeCloud((0,0,0),(10,10,20),rMean=0.5)

# add the sphere pack to the simulation
sp.toSimulation(material=Mat2)

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
InteractionLoop(
# handle sphere+sphere and facet+sphere collisions
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.0,label='integrator'),

# call the checkUnbalanced function (defined below) every 2 seconds
PyRunner(command='deposition()',realPeriod=2,label='checker'),
# call the addPlotData function every 200 steps
]

O.dt=.5*PWaveTimeStep()

def deposition():
if O.iter<100: return
if unbalancedForce()>0.1: return
O.bodies.append(utils.wall(position=(0,0,20),axis=2,sense=0,material=Mat1))
global plate
plate=O.bodies[-1]
plate.state.vel=(0,0,-vel)
checker.command='stopPlate()'

def stopPlate():
if plate.state.pos[2]<15:
plate.state.vel*=0
integrator.gravity=(0,0,+9.81)
checker.command='stop1()'

def stop1():
if unbalancedForce()<.02:
O.pause()

)

plot.plots={'i':('UnbalancedForce',None,'Porosity','Density',
)}

plot.plot()

O.saveTmp()

Thank you guys!! =)

 Jan Stránský (honzik) said on 2020-02-06: #6

> However I cannot make the top stop at a desired height.
> def stopPlate():
> if plate.state.pos[2]<15:
> plate.state.vel*=0
> integrator.gravity=(0,0,+9.81)
> checker.command='stop1()'

here the change of checker.command is done in the very first run (and the plate is then not stpped), maybe put it inside the if condition?

cheers
Jan

 Paolo (p4olo) said on 2020-02-06: #7

Oh, the indentation.... a noob's mistake.

It works very good now!

I have just a final question and I'm sorry for the off topic. Could you please tell me if it is possible to define the order of magnitude/scale of the axes in the plot?
i.e. I want to plot the density/porosity vs iteration like 0-100 (%) but the y axis is in automatically set in scale 1e1 so it is 0-10 (-)

Tank you all so much for the help and the time!!

 Jan Stránský (honzik) said on 2020-02-09: #8

Next time please open a new question for a new topic [1].

if the problem is the units, than simply save the value in % (100 x higher).
If it is about setting specific range, I don's know (then please follow the first line of this post)

cheers
Jan

 Launchpad Janitor (janitor) said on 2020-02-22: #9

This question was expired because it remained in the 'Open' state without activity for the last 15 days.