How to reset the repulsive forces at the beginning of a compression test

Hi Damien,

it seems like a bug in RegularHexa function. After adding

O.step()
for i in O.interactions: print i.geom.penetrationDepth

at the end of your code, I got values like 1.73472347598e-18, which is ok
and is the result of rounding error, but also values 0.000757886244659,
which is definitely not ok. Feel free to report the bug.
Cheers
Jan

2015-05-12 17:16 GMT+02:00 GERARD Damien <
<email address hidden>>:

> New question #266828 on Yade:
> https://answers.launchpad.net/yade/+question/266828
>
> Hello,
> This question is related to question #264456 :
> https://answers.launchpad.net/yade/+question/264456
> I created a spherical assembly made up of smaller spheres bounded with
> cohesion which I submit to compression. The ojective of the simulations is
> to reproduce failure in the spherical assembly induced by successive
> decohesion.
> Due to initial interpenetration, there is some repulsive forces that broke
> the cohesion at the beginnig of the test (If I increase the gap between
> balls there are not enough contacts between spheres.
>
> So I though to reset the repulsive forces to solve this problem, by using
> unp equal to penetration depth.
> But I couldn't manage to succeed in this task.
>
> Is anyone can help me to do this?
>
> thanks.
>
> Regards, Damien
>
> this is piece of my code (I left only the creation of the sphere
> assembly-the problem appears even without loading this assembly):
>
> from yade import pack
> from math import *
>
> cohN=50e6
> cohT=35e6
> E=50e9
> nu=0.3
> d=2310
> phi=35
> rayon_s=0.005
> rayon_g=0.05
>
>
> O.materials.append(CohFrictMat(young=E,poisson=nu,density=d,frictionAngle=radians(phi),normalCohesion=cohN,shearCohesion=cohT,label='sol'))
>
>
> sp=pack.regularHexa(yade._packPredicates.inSphere(center=(0,0,0),radius=rayon_g),gap=0,radius=rayon_s,material='sol')
> O.bodies.append(sp)
>
> 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()),
> NewtonIntegrator(damping=0.3)]
>
> O.dt=utils.PWaveTimeStep()
> O.usesTimeStepper=True
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hi,
For arbitrary initial positions you can define the the equilibrium displacement via i.phys.unp like this:
i.phys.unp = i.geom.penetrationDepth #I didn't check the syntax but it should be nearly this
Then initial positions will correspond to zero forces.

Should be fixed:

https://github.com/yade/trunk/commit/d384991b5866aaa3a608a54e592e0feca7d8d225

Note that you may face some issues using the assignment phys.unp = geom.penetrationDepth, because of an "offset" (during the computation loop) in the computations / use of these two variables.
See the example below, where some movements (vanishing here after some oscillations) appear.

If, in your case, such oscillations would be annoying, e.g. having an amplitude big enough to break the cohesive bonds, you may just perform one cycle with any movement frozen (dynamic = False).
So that the initial offset between unp and penetrationDepth can be cancelled (uncomment the two corresponding lines in the example, and you won't get any movement)

O.materials.append(CohFrictMat(young=50e9,poisson=0.3,frictionAngle=radians(35),normalCohesion=1e10,shearCohesion=35e10))

O.bodies.append(sphere(center=Vector3(0,0,0),radius = 1,fixed=1))
O.bodies.append(sphere(center=Vector3(1.5,0,0),radius = 1))

O.engines=[
    ForceResetter()
    ,InsertionSortCollider([Bo1_Sphere_Aabb()])
    ,InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D()],
        [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False)],
        [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()])
    ,NewtonIntegrator()
    ,PyRunner(command='saveData()',iterPeriod = 1)
    ]

O.dt = 1e-5
from yade import plot

def saveData():
    plot.addData(uSph = O.bodies[1].state.pos[0]
                 ,fN = O.forces.f(1)[0]
                 ,it = O.iter
                 )

plot.plots={'it':'fN',' it':'uSph'}
plot.plot()
yade.qt.View()

#O.bodies[1].dynamic=False
O.step()
O.interactions[0,1].phys.unp = O.interactions[0,1].geom.penetrationDepth
#O.bodies[1].dynamic=True
O.run(1000,1)

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

Hi Jérôme
I tried the script with Yade 1.12 but below error was showed in terminal:

Running script checkRepulsiveForcecheckRepulsiveForce_266828_4.py
Traceback (most recent call last):
  File "/usr/bin/yadedaily", line 180, in runScript
    execfile(script,globals())
  File "checkRepulsiveForcecheckRepulsiveForce_266828_4.py", line 26, in <module>
    O.interactions[0,1].phys.unp = O.interactions[0,1].geom.penetrationDepth
IndexError: No such interaction
[[ ^L clears screen, ^U kills line. F12 controller, F11 3d view (use h-key for showing help), F10 both, F9 generator, F8 plot. ]]

As you may have figured out, the error says that there is no interaction between the 2 spheres. Thus, it is not possible to access either the phys or the geom of this not-existing interaction, then the crash of the script.. This is for the description of the problem.

Unfortunately, I can not imagine what the reason could be. Your problem here is that Yade does not get there should be an interaction between two overlapping spheres...

* did you run the exact same example that I wrote above ? Copying-pasting it in a text file (with "O.materials" on line 1), I find the line 26 is the "plot.plots" line, and not the "O.interactions[0,1]...." one
* do you confirm the error appears at first iteration ? (what is the result typing "O.iter" ?)

Hi Jerome,

I have been reading your comment on the introduction of an equilibrium distance using cohFrictMat.
I am not sure I understood properly what you point out : my interrogation comes from reading this previous thread
 --> https://answers.launchpad.net/yade/+question/295035,

where the problem of initial offset or undesired movements was not raised.

I do not understand why there should be a problem in the first computational loop, is it not possible to loop over all interactions in order to set i.phys.unp = i.geom.penetrationDepht BEFORE entering the actual computation loop ? Or does Yade need to actually perform a first loop in order to set initial values of the simulation ?

Plus, I was also reading the thread https://answers.launchpad.net/yade/+question/269724

 discussing whether Yade would match the BPM model features, and you mentioned that in JCFpmMat, the initial distance between particles is automatically set as the equilibrium distance ? So there is no need to use the initD function in order to achieve a zero stress state for initial overlapping particles using that contact law ?

Regards

Hi,
>"I do not understand why there should be a problem in the first computational loop, is it not possible to loop over all interactions in order to set i.phys.unp = i.geom.penetrationDepht BEFORE entering the actual computation loop ?"

Because before first iteration bodies are there but no interactions yet (there are exceptions to this rule, maybe you have one in mind).
There are many ways to solve this problem, namely.
1/ set O.dt=0 and run one iteration in order to create the interactions in a fixed system, then set unp.
2/ execute "O.engines[1](); O.engines[2]()" (collision detection + interaction loop), then set unp
3/ wait until Janek get the feature implemented ;)
Bruno

Bruno already answered the first part of your question.
If I want to rephrase, I can say that most often the 1st computation loop is required to create the interaction network, thus no, it is not possible to set i.phys.unp = i.geom.penetrationDepth before this computation loop..
See Bruno's answer for ways to deal with this feature, or also another (less convenient ?) way in [1]

For the 2d part of your question and JCFpm, I confirm this initD attribute is automatically correctly defined at the interaction creation, see [2].
So, during a JCFpm simulation, you would actually use this initD attribute, but without knowing it.

Thank you for browsing the archives ! :-)

[1] https://yade-dem.org/doc/user.html#individual-interactions-on-demand
[2] https://github.com/yade/trunk/blob/e4e757f2e98a620e3177b7a36a1d10f69f6a6a28/pkg/dem/JointedCohesiveFrictionalPM.cpp#L39

@Jérôme: [1] is helpless, how would you use it?

Following #6 in https://answers.launchpad.net/yade/+question/233232... ;-)

It was for creating interactions between distant objects, _provided that_ the list of interactions to be created is known in the first place. I don't see how to apply this in the present case (maybe it was also irrelevant in the context of #6? not sure).

Well,
Thanks Bruno and Jerome, I understood this issue. I am still working on understanding how Yade fully works.

that's why I read archives !

Regards,

I pick option 3), should be ready in few weeks :)