Simple compression on an sphere assembly

Hi,

Your description still sounds much like a numerical instability problem. To be sure it is not time step related, a first (or second since you tried already to reduce dt) attempt could be to introduce a GlobalStiffnessTimeStepper engine in your O.engines. It is very efficient to set a convenient time step

Hello,
I don't see a real problem in the script you sent. Even multiplying both
timestep and velocity by 10 gives a steady compression of the cloud
(You can also decrease damping to a more reasonable value, 0.7 is huge).
Where do you see a problem?

Bruno

Hello,
Firstly thank you for the advice
Even adding the option timestepper, the beginning of my test ball explodes,
I like to know that it is the problem that the ball explodes when the gap is almost zero.
Thank you in advance

I can't reproduce any explosion with the provided script (with yade 1.12.0-94) . I don't understand the problem, sorry.

Hello,
I would like to make additions compared to my old message.
By testing high Young's modules, the balls go in all directions. Here is my script.
What should I do?
I thought it was necessary to reset to zero contact forces at the first time step but I do not know what commands I have to write in my script. Do you think it's the right solution? and what are these commands?
Moreover, what are the commands to know the value of the contact forces for each interaction / contact?
Thank you in advance.
--------------------------------------------------------------
from yade import pack
from math import *
from yade import plot
from yade import qt

vitesse=0.1
cohN=50e5
cohT=35e5
E=50e9
nu=0.3
d=2310
phi=35
rayon_s=0.005
rayon_g=0.05
rayon_boite=rayon_g+0.0001
espace=0.0
gravity=0

O.materials.append(CohFrictMat(young=E,poisson=nu,density=d,frictionAngle=radians(phi),normalCohesion=cohN,shearCohesion=cohT,label='sol'))
O.materials.append(CohFrictMat(young=10e13,poisson=0,density=10000,frictionAngle=radians(10),normalCohesion=0,shearCohesion=0,label='walls'))

mn,mx=(-rayon_boite*2,-rayon_boite,-rayon_boite*2),(rayon_boite*2,rayon_boite,rayon_boite*2)
walls=utils.aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.regularHexa(yade._packPredicates.inSphere(center=(0,0,0),radius=rayon_g),gap=espace,radius=rayon_s,material='sol')
O.bodies.append(sp)

for i in O.interactions:
 i.phys.unp = i.phys.penetrationDepth

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1),Bo1_Box_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=1),Ig2_Box_Sphere_ScGeom6D()],
        [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False),Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
        ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=5,timestepSafetyCoefficient=0.8, defaultDt=utils.PWaveTimeStep()),
    TranslationEngine(
    translationAxis=(0,-1,0),
    ids=[3],
    label='tl',
    velocity=vitesse
    ),
    NewtonIntegrator(damping=0.3, gravity=[0, gravity, 0]),
]

O.dt=utils.PWaveTimeStep()
O.usesTimeStepper=True

qt.Controller()
qt.View()

O.saveTmp()

Hello Damien,

Indeed, in this new version of script, I understand you re unhappy with your spherical aggregate imploding. Which was not the case in the first version of the script, where, like Bruno, I do not see any problem.
Do not hesitate to state such things, so that we understand more what's going on along this thread.

If one script is "OK", and the other is not, the first thing to solve your problem is to check what differences exist in both scripts ! ;-)

You're telling us about higher "Young" parameters. It is not the reason.
Just put in the second script E=50e7 (like in the first one), and you will still get this annoying behaviour.
Note that with numerical instabilities problem, spheres generally fly to the moon, while here they stay in the range of the graphical view. So, it is probably (surely in fact) not an instability problem related to excessively high stiffnesses. Which is logical since you're using this GlobalStiffnessTimeStepper which "is very efficient to set a convenient time step", as I said before.

On the other hand, you do not tell us about other differences between both scripts
- you changed the "gap" value of pack.regularHexa
- you added these lines
 for i in O.interactions:
 i.phys.unp = i.phys.penetrationDepth

Why ? And why did you want to reset zero contact forces (in fact I have an idea, but it would help that you state it yourself ;-) )

For your other question (better open different Launchpad questions when you have different questions), the contact force of a given interaction "i" is accessed with
i.phys.normalForce. "i" being a member of O.interactions (see https://yade-dem.org/doc/yade.wrapper.html#interactioncontainer)
Or, more simply, through the "Interaction" tab of "Inspect" window of Yade graphical interface.

Jerome

PS : if you are yourself lost to identify changes between scripts, you may use "kompare" a nice software
PPS : providing scripts hoping that others will test and dig into them is a great favour you are asking. You will have more chance to get this favour removing as much stuff as possible from your script.
E.g. the definition of a zero gravity, and of a Translation Engine, which is not related to your problem

Hello Jérôme,
Thank you for your answer and your time.

I explain you my different changes :
- It is true that I changed the "gap" value of pack.regularHexa to zero because as I tried to evaluate the influence of this parameter on the results, to find an issue of my problem. It is a parameter which can be tuned on.

-The lines " for i in O.interactions: i.phys.unp = i.phys.penetrationDep" have been added because I though that reset zero contact forces at the beginning could solve my problem, as I read in other answers on this forum (https://answers.launchpad.net/yade/+question/257058).
But I couldn't manage to succeed in this task.

If you need more clarifications, do not hesitate to tell it.
Regards, Damien

Thanks Jérôme Duriez, that solved my question.