Using ForceEngine to break a single JCFpm bond

Asked by Yaniv Fogel on 2019-09-09

Hello,

I'm a new user to YADE and I'm in the process of learning about the JCFpm model. I'm currently trying to break a bond between two particles using the ForceEngine. As described in the code bit below, I'm trying to apply and increasing force until a break occurs between the two particles. The bottom particle is fixed.

If I got the documentation right, the force when the breakage should happen is:
FnMax = tensileStrength * pi*Rmin^2 [1],[2] -> FnMax = 1e6 * np.pi * pow(0.005,2) = 78.539 N

First question -
For the first trials, I tried using addF manually, using the code below (after running the script). I found that the bond breaks around O.forces.addF(1,(0,0,700)). (I first assumed that addF uses Newtons, but now I'm not so sure...). So, did I calculate FnMax correctly or not?

Second question -
After that, I tried using the ForceEngine (same code, just uncomment #ForceEngine & #PyRunner), increasing the force on the particle, but I can't get it to break. Am I using the engine correctly?

Third question -
I also looked at the InterpolatingDirectedForceEngine [3], where there is a magnitude for the force, but what is the difference between force and magnitudes in the engine (Again, I tried playing with it, with increasing magnitudes, but I still couldn't get it to break)

I would like your help. Sorry if this is too long or if my understanding of vector physics is a bit lacking :)

### Forces Test 1, YADE ###

from yade import pack, plot
import numpy as np

yade.qt.Controller()

# Define geometry
r = 0.005
intR = 1.0

# Define material
idRockTest = O.materials.append(JCFpmMat(type=1,
                                         young=30e9,
                                         density=2500.0,
                                         poisson=0.1,
                                         frictionAngle=np.radians(18.0),
                                         tensileStrength=1e6,
                                         cohesion=1e6,
                                         label='Rock'))

# add spheres
O.bodies.append([
    sphere((0,0,0),r,material='Rock',color=(0.019, 0.529, 1),fixed=True),
    sphere((0,0,2*r),r,material='Rock',color=(1, 0, 0))
    ])

# Simulation loop
O.engines=[
 ForceResetter(),

    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),

    InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom')],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(label='interactionLaw')]
 ),

    # Apply gravity force to particles.
 NewtonIntegrator(
        gravity=(0,0,-9.81),
        damping=0.2),

    # Apply force on top particle
    #ForceEngine(ids=[1],force=Vector3(0,0,0),label='fEngine'),

    #PyRunner(command='fEngine.force[2]+=100',realPeriod=1,initRun=True),

    PyRunner(command='recorder()',realPeriod=1,initRun=True),
]

# Define time step
#O.dt = 0.5*utils.PWaveTimeStep()

def recorder():
    plot.addData({
        'i':O.iter,
        'F':O.forces.f(1).norm(),
        'T':O.forces.t(1).norm()})

plot.plots={'i':('F','T')}
plot.plot()

O.dt=0.

#### manage interaction detection factor during the first timestep and then set default interaction range
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

# Define time step
O.dt = 0.5*utils.PWaveTimeStep()

O.run()

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.JCFpmMat.tensileStrength
[2] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.JCFpmPhys.crossSection
[3] https://yade-dem.org/doc/yade.wrapper.html?highlight=interpolatingdirectedforceengine#yade.wrapper.InterpolatingDirectedForceEngine

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2019-09-25
Last query:
2019-09-25
Last reply:
2019-09-25
Jérôme Duriez (jduriez) said : #1

Hi,

If your goal is to learn about JCFpm (and not Force Engine in itself), I'd advice to better perform a displacement controlled test, where one particle is pulled from with the other with a constant velocity.
(both particles being non dynamic, thus)

And measure the interaction force as long as the interaction lasts.

Might be easier to interpret.

Regarding FnMax, see also Eq. (2) in https://www.sciencedirect.com/science/article/pii/S0013794415007225 where t therein is YADE's JCFpmMat.tensileStrength

Yaniv Fogel (yanivfgl) said : #2

Hi Jérôme,

Thanks for the tips, but I am interested on the ForceEngine, rather than the JCFpm.
I will do the displacement test just to verify the values. But again, I try to do it using a force... and I'm still unable to solve this...

I imagine the ForceEngine as "tying" a string to the top particle, and applying tension gradually... As I see it in my imagination, at certain tension, it should break, and the partilce should start moving... but that doesn't happen...

Is the way I "imagine" this action is wrong?

Yaniv

Jérôme Duriez (jduriez) said : #3

As for ForceEngine, it's maybe better to imagine it as injecting, in the sum of forces acting on bodies "ids", a new force.

What happens next comes from the dynamics of the system.

In particular, you may relate the ForceEngine force with the interaction force only in special cases (only two forces on body named 1 here : ForceEngine.force, and the interaction force ; combined with quasi static conditions where these two forces are just opposite).

PS : your script is too long for me to give a real look at it :-)

Jan Stránský (honzik) said : #4

Hello,

> First question -
> For the first trials, I tried using addF manually, using the code below (after running the script). I found that the bond breaks around O.forces.addF(1,(0,0,700)). (I first assumed that addF uses Newtons, but now I'm not so sure...). So, did I calculate FnMax correctly or not?

as suggested by Jerome, there is some dynamics in the system. You area inreasing the applied force and the particle accelerates, but the force in the interaction would be (most likely) less than the force applied (determined by displacement).
Then, naturally, the break is "delayed" from the force applied. So apart from applied force, you can also track the data from interaction itself.

> Second question -
> After that, I tried using the ForceEngine (same code, just uncomment #ForceEngine & #PyRunner), increasing the force on the particle, but I can't get it to break. Am I using the engine correctly?

The ForceEngine and PyRunner is commented (?) and therefore has no effect..

cheers
Jan

Yaniv Fogel (yanivfgl) said : #5

Hi Jan and Jérôme,

Thanks for the remarks. I understand what you mean about the dynamic state of the system.

I tried investigating how the forceEngine and interaction forces are related. I would expect that as the forceEngine force increases, the interaction force will decrease (as they are opposite), but like before, there is no reaction between the forceEngine and interaction forces between the particles...

As you will see while you run the script provided below ( I tried to see how I can minimize it anymore, but I think this is the most minimal way...), the force acting on top particle (F1) is increasing gradually, while the force acting on the bottom particle (F0) is the gravity force, and is logically equal to the interaction force between the particles. The engineForce has no effect on the system... :/

(There are some "steps" sometimes in the plot of F1 between interactions 0 and 0.2e7, when the force reaches around 70-120, which does correlates to a certain degree with the FnMax value I calculated ... but again, other than that, there is no effect...)

#### Updated forceEngine Code ####

 from yade import pack, plot
import numpy as np

yade.qt.Controller()

# Define geometry
r = 0.005
intR = 1.0

# Define material
idRockTest = O.materials.append(JCFpmMat(type=1,
                                         young=30e9,
                                         density=2500.0,
                                         poisson=0.1,
                                         frictionAngle=np.radians(18.0),
                                         tensileStrength=1e6,
                                         cohesion=1e6,
                                         label='Rock'))

# add spheres
O.bodies.append([
    sphere((0,0,0),r,material='Rock',color=(0.019, 0.529, 1),fixed=True),
    sphere((0,2*r,0),r,material='Rock',color=(1, 0, 0))
    ])

# Simulation loop
O.engines=[
 ForceResetter(),

    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),

    InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom')],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(label='interactionLaw')]
 ),

    # Apply gravity force to particles.
    NewtonIntegrator(
        gravity=(0,-9.81,0),
        damping=0.2),

    # Apply force on top particle
    ForceEngine(ids=[1],force=Vector3(0,0,0),label='fEngine'),

    PyRunner(command='fEngine.force[1]+=10',realPeriod=1,initRun=True), # apply gradual increasing force on top particle

    PyRunner(command='recorder()',realPeriod=1,initRun=True),
]

# plot forces and position
def recorder():
    plot.addData({
        'i':O.iter,
        'F1':O.forces.f(1).norm(), # force acting on top particle
        'Pos[z]':O.bodies[1].state.pos[1], # position (y) of top particle
        'F0':O.forces.f(0).norm(), # force acting on bottom particle
        'nInter[0,1]':O.interactions[0,1].phys.normalForce[1]}) # normal interaction force between particles

plot.plots={'i':('F1'),'i ':('F0'),'i ':('Pos[z]'),'i ':('nInter[0,1]')}
plot.plot()

O.dt=0.

#### manage interaction detection factor during the first timestep and then set default interaction range
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

# Define time step
O.dt = 0.5*utils.PWaveTimeStep()

O.run()

Jérôme Duriez (jduriez) said : #6

If you're interested into ForceEngine, I would, if I were you, first of all investigate ForceEngine alone before combining ForceEngine and JCFpm material.

A minimal working test to study ForceEngine's behavior and check you're correctly using it could be this one : (Python 2 in the print line, sorry...)

########### Beginning #####
# one sphere with weight imposed by ForceEngine for sake of ForceEngine investigation (use NewtonIntegrator.gravity otherwise):
O.bodies.append(sphere((0,0,0),1))
O.engines = O.engines[0:-1] + [ForceEngine(ids = [0],force = Vector3(0,0,-O.bodies[0].state.mass*10))] + [O.engines[-1]] # introducing ForceEngine to apply weight, taking advantage of default O.engines

O.engines[5].damping = 0 # get rid of this numerical damping trick for a correct interpretation of the results / theory

O.dt = 0.01
O.dynDt = False
O.run(100,True)

print '\nAfter 100 iterations id est 1 second of imposed free fall, the body is at z =',O.bodies[0].state.pos[2],'versus -5 expected\n'
########### End ########

You should get good results here / get more confident with ForceEngine.

Best Jan Stránský (honzik) said : #7

Hi,

a few minor notes:

> PyRunner(command='fEngine.force[1]+=10',realPeriod=1,initRun=True),
> PyRunner(command='recorder()',realPeriod=1,initRun=True),

use iterPeriod or virtPeriod

> O.run()

use something like O.run(1000,True)

> The engineForce has no effect on the system... :/

the most important thing: put ForceEngine before NewtonIntegrator

cheers
Jan

[1] https://yade-dev.gitlab.io/trunk/user.html#base-engines, "if it apply force, it should come before NewtonIntegrator ..."

Yaniv Fogel (yanivfgl) said : #8

Thanks Jan Stránský, that solved my question.

Yaniv Fogel (yanivfgl) said : #9

Hi all,

Jérôme, Thanks for the tips. I will dig deeper into the ForceEngine.

Jan, That last important tip solved the problem :) I can't believe I missed that line in the documentations.
I can say that the JCFpm bond broke at the calculated FnMax :)