GSTS gives an unstable dt?

Asked by Alessandro Finazzi on 2019-03-08

Hello everyone,

I'm Alessandro and I'm a new Yade user. I have to simulate the behaviour of a given set of spheres (radii, material, accelerations are given) inside a cylinder. I use yadedaily, xenial version.

The given set of accelerations translates into amplitudes of the HarmonicMotionEngine up to 7 mm. And here I find my problem: up to 1-2 mm all is ok, no unstability problems. But when I try, for example, with an amplitude of 3.5 mm, the spheres begin to overpass the facets of the cylinder and, as far as I know, it could be an issue linked to instabilities.
Firstly I tried to use a smaller dt, then I read about GSTS and decided to move to that. It wolud be perfect for my purpose: in the first part of the script I need a sort of "gravity deposition" and so a bigger dt would be nice, then the GSTS should assure the right dt in order to avoid instabilities (if I understood correctly).

But using GSTS, even with a very small safety coefficient, the spheres keep overpass the facets. I have to admit that with a lower coefficient the simulation seems to be improved, but the problem is not solved.

So: do you think that's an instability issue? Is the decreasing of the safety coefficient the only way to improve it? And shouldn't GSTS give a pretty good dt even without this very small coefficient?

I've used both CoundallStack (now commented) and ViscElPhys, and my problem is present in both cases.

Here I attach my script:

### MOVIMENTO ARMONICO ALTO BASSO

from yade import pack,utils,plot,ymport,export,qt,timing

## Physical Parameters
density=7950
en = 0.5
es = 0.5
cyl_r = 12.7
cyl_h = 24.6
sp_p_r = 1.5
fill_size = 11.5
young = 200e9
poisson = 0.305

## Creation Cylinder
facetMat=O.materials.append(ViscElMat(density=density,young=young, poisson=poisson,en=en,et=es))
Cylinder=O.bodies.append(geom.facetCylinder((0,0,cyl_h/2),cyl_r,height=cyl_h,segmentsNumber=15,wallMask=7,material=facetMat))

## Creation Spheres
sphereMat=O.materials.append(ViscElMat(density=density,young=young,poisson=poisson,en=en,et=es))
sp=pack.SpherePack()

sp.makeCloud((-(cyl_r-2),-(cyl_r-2),1.5),(cyl_r-2,cyl_r-2,fill_size-1.5),rMean=sp_p_r,num=50,distributeMass=False)
sp.toSimulation(material=sphereMat)

sp2=pack.SpherePack()
sp2.makeCloud((-(cyl_r-2),-(cyl_r-2),fill_size),(cyl_r-2,cyl_r-2,cyl_h),rMean=3,num=3,distributeMass=False)
sp2.toSimulation(material=sphereMat)

## Engine
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
      [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
      #[Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_ViscElPhys_Basic()]
      #[Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   GlobalStiffnessTimeStepper(timeStepUpdateInterval=1,label='tss',timestepSafetyCoefficient=0.001,viscEl=True),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.3),
   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker')
   ]

def checkUnbalanced():
 if O.iter<15000: return
        if kineticEnergy()>18000000: return
        O.engines=O.engines+[HarmonicMotionEngine(A=[0,0,3.5],f=[0,0,60.0],fi=[0,0,pi],ids=Cylinder), PyRunner(command='addPlotData()',iterPeriod=10)]
 checker.command='nothing()'

def nothing():
        return

def addPlotData():
 Ek,maxId=kineticEnergy(findMaxId=True)
 plot.addData(i=O.iter,t=O.time,Ninter=O.interactions.countReal(),total=O.energy.total(),maxId=maxId,**O.energy)

O.trackEnergy=True
O.saveTmp()

plot.plots={'i':('Ninter',),'t':(O.energy.keys,None,'total')}

plot.plot()

qt.Controller()
qt.View()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Bruno Chareyre
Solved:
2019-03-08
Last query:
2019-03-08
Last reply:
2019-03-08
Jérôme Duriez (jduriez) said : #1

Hi,

Which bodies have their displacements/positions controlled ?
the spheres (as your question says in the beginning) or some boundary facets (HarmonicMotionEngine seems to apply to those only) ?

Note that displacement controlled bodies usually should have state.dynamic = False. I also think it may be simpler to directly modify their state.vel, instead of using (old) Engines whose exact behavior may be more complicated to predict for a new user. (sorry about that..)

Generally speaking, my first guess would be:

No, it's not an instability issue (because GSTS seems to me to do its job very well).
Maybe your spheres and facets move so fast that they can cross each other in one timestep, without noticing it. This would be also time step-related obviously, but not in a divergence sense.

Unstable dt give weird results in the sense they go to infinity (= divergence)

Fabio (fabiogabri) said : #2

Hi,
I don't check your code but I guess in some (rare) cases it could be a problem related to the number of spheres you used, the bulk porosity, and their size in comparison to the amplitude of HarmonicMotionEngine. Let suppose to use two particles vertically aligned, and a plane at the base. The upper one is displacement controlled. If you move the upper one with an amplitude exceeding 2*radius the particle at the bottom could pass through the plate (or the upper sphere). The time step can be small but I suppose it doesn't matter.

One more option: it could be that timestep is ok wrt both numerical stability and velocity vs. diameter, with particles going through the plate due to the magnitude of forces.
young = 200e9 sounds stiff but you are speaking of moving 25m of steel with +/- 3.5m amplitude at a frequency of 60Hz, which is beyond current technology .
Check your units.

Funny question: you describe "A=[0,0,3.5]" as an amplitude of 3.5mm. How would you change the script for an amplitude of 3.5m?
What about "density=7950"? micropounds per cubic foot?

Hi Jérôme and Fabio,
thanks for your answers.

Just to answer to Jérôme question and to explain in a more detailled way the simulation that I need: I have a sort of mixer, filled with spheres, that is shaken up and down at very high accelerations (from 50 g to 100 g, as peak). And I have to study energies, number of interactions and so on of the system.
In short: just the facets composing the cylinder will have their position controlled.

From your answers I get that isn't an instability issue, and so GSTS isn't resposible for this behaviour.
I'll try to modify the script and use state.vel instead of the HarmonicMotionEngine, let's see if it solves the problem.

If you have any other suggestion to avoid this "crossing" that you described, let me know.

Hi Bruno,
you totally got the point.
I was messing around with very small time steps while I had taken m for mm. I corrected the script and now it seems to work perfectly. Beyond any doubt the problem was given by the magnitude of forces.

Thanks to Jérôme and Fabio for your explanations and your fast answers, and thanks another time to Bruno that solved the problem pointing out my mistake.

Thanks Bruno Chareyre, that solved my question.

Fabio (fabiogabri) said : #8

Your initial cloud seems to be external to the cylinder...