Yade

Unbalanced simulation

Asked by Nolan Dyck

Hello again,

I am having some trouble achieving a stable packing using small particle sizes (isotropic compression, strain controlled). Using the script below I can get stable packings when the particle sizes are big enough (e.g. set SC=1). But when I decrease the size of the particles (e.g. set SC=1e-6) the simulation starts off fine but then the particles begin flying around and the simulation never reaches a stable packing. Now I have tried playing around with several parameters to get it to work: time step, material density, particle elasticitiy, verletDist, damping, maxUnbalanced, globUpdate, maxStrainRate. Decreasing the timestep seems to have a small effect on stability, but it yeilds unbearable simulation times. Note that the number of particles in my simulation is relatively small (I will control it later) and I don't want to decrease it anymore. What do I need to do to make stable packings for small particle sizes?

As a side note, when I decrease the particle sizes I will get an error telling me that some particles are bigger than half the cell period, even though its perfectly clear from the OpenGL viewer that this is not the case (is this a bug?). For now I have set the allowBiggerThanPeriod flag to true.

I have also looked at this post:
https://answers.launchpad.net/yade/+question/228886
but I don't really understand how the kinetic energy problem was solved.

Thanks,
Nolan

CODE:

from yade import pack,plot,export

r=2.12e-1
p=0.70
R=2.63155
t=0.065

print 'r = ',r
print 'R = ',R

SC=1e-6
print 'SC = ',SC

O.dt=5e-7
print 'O.dt = ', O.dt

O.materials.append(FrictMat(young=1e9,poisson=.35,frictionAngle=radians(0),density=1e9,label='spheres'))
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(R*SC,R*SC,R*SC),rMean=r*SC,rRelFuzz=0,periodic=True)
sphVol=sp.relDensity()*(R*SC)**3
tStrain=-1
StrainR=0.1

sp.toSimulation()
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=0,allowBiggerThanPeriod=True),
   InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_ScGeom_FrictPhys_CundallStrack()]),
   NewtonIntegrator(damping=0.95),
   PeriTriaxController(goal=(tStrain,tStrain,tStrain),stressMask=0b000,maxUnbalanced=1e-10,globUpdate=5, doneHook='reachedTargetStrain()',label='triax',maxStrainRate=(StrainR,StrainR,StrainR),dynCell=True),
   PyRunner(command='checkPorosity(REVl,sphVol)',iterPeriod=50)
]
O.run()

def reachedTargetStrain():
   print 'Reached Target Strain'
   O.pause()

def checkPorosity(REVl,sphVol):
 REVl = O.cell.size[0]
 REVvol = REVl**3
 por = sphVol/REVvol
 #print 'porosity = ',por
 if(por > p):
  print 'Target porosity reached'
  print 'Finished Simulation'
  print 'porosity = ',por
  print 'REVl = ', REVl
  O.pause()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Bruno Chareyre
Solved:
Last query:
Last reply:
Revision history for this message
Hien Nguyen (giahien) said : 		#1

Hi, just add the comment because you quoted my previous question :)
In fact in the past my model have very high value of kinetic energy because I fill a very small particles in a big box (after that I make the walls move to compress).
I try to reduce the box's size (the utils.aabbWalls() functor) and modify the initial particles' size (the rMean and porosity in O.makeCloud) and each time modifying, I run my script again until I have the acceptable reponse of the kinetic energy.
I don't know is there any other reason, but for my problem before, the initial size of the box, the cloud of particles and the maxMultiplier+finalMaxMultiplier (if you choose the growing sphere method) or the speed you choose for the wall (max_vel, if you use moving wall method) are very important parameters, try to modify them until you have a nicer kinetic energy.
I will follow this question :)

Revision history for this message
Best Bruno Chareyre (bruno-chareyre) said : 		#2

Nolan,
Explicit time integration is conditionally stable. Smaller particles require smaller timesteps.
I suggest using GlobalStiffnessTimeStepper to set dt automatically, as in the example triaxial scripts.

https://www.yade-dem.org/doc/formulation.html#stability-considerations

Revision history for this message
Nolan Dyck (ndyck) said : 		#3

Thanks Bruno Chareyre, that solved my question.

Revision history for this message
Nolan Dyck (ndyck) said : 		#4

Nguyen,

Thanks for clarification about that thread. Unfortunately I realized that we are using different engines so I am not controlling the same parameters.

Bruno,

Thanks for the advice. I'm running a test simulation and the GlobalStiffnessTimeStepper seems to be working well (although it is slow). I'm just happy that it seems to promise a successful simulation rather than an unstable one.

That being said, I have another issue. As you know I am working on bubble interactions. Right now those interactions have no stiffness value so I will have to come up with a substitute. It seems a little complicated to reformulate the GlobalStiffnessTimeStepper but I could linearize the bubble constitutive law and then have it report an effective contact stiffness to the time stepper. I can write the code to calculate the effective stiffness easily enough, I just need to figure out how to report it to the GlobalStiffnessTimeStepper. Also, is this a good approach? Do you have any better (or simpler) ideas?

Thanks,
Nolan

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#5

>GlobalStiffnessTimeStepper seems to be working well (although it is slow)

GSTS is slow or your simulation is slow?
When dt is set by GSTS, you can always pause the simulation and set dt=2*dt manually (in controller window) to see if GSTS is choosing an excessively small dt. In my experience it never happens, and dt=2*dt should explode immediatly (unless you changed the default safety factor to something smaller than 0.8).

Classical tricks for speedup: increase mass or decrease stiffness.

>I could linearize the bubble constitutive law and then have it report an effective contact stiffness to the time stepper.

Exactly. All you need is to reflect the current (or expected maximum, if you don't want to update all the time) linearized stiffness in the kn/ks variables of the contacts. GSTS will use them to give a timestep.

Revision history for this message
Nolan Dyck (ndyck) said : 		#6

I was just acknowledging that the simulation was slow (compared to large particles).

Thanks for your help again, Bruno.