Apply cohesive interaction on JCFpmMat

Hi,

Yes, there is an easier way. Juste define the material parameters in such lines:

def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson = ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
def wallMat(): return JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))

Second, make sure that you have these lines in your engines list:

...
InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
),
...

with intR a value that you previously determined to ensure a coordination number K equal to the value you want for the packing you use (K=10 for sandstone).

Then, the key parameter is: cohesiveTresholdIteration which defines the iteration up to which you want to create bonds between particles (the first one here).

Thus, just by runnning

O.step()

the bonds will be created.

After, you may want to reset the interaction radius to 1 by typing

SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

You can check that the coordination number K corresponds to the one you want by:

numSSlinks=0
numCohesivelinks=0
for i in O.interactions:
    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
      numSSlinks+=1
    if i.phys.isCohesive :
      numCohesivelinks+=1
print "| K=", 2.0*numCohesivelinks/numSpheres

and then run your simulation (triaxial loading).

Luc

Thanks for Luc's solution. But I still have a problem:
In triaxial test, the confining pressure is applied through internal compaction. And the isotropic state is reached at a certain step(10000 for example).
The method provided by Luc set enlargeFactor and cohesion at the first step. Then how should these two parameters be set at 10001?

--

At 2017-03-13 20:52:53, "Luc Scholtès" <email address hidden> wrote:
>Your question #556907 on Yade changed:
>https://answers.launchpad.net/yade/+question/556907
>
> Status: Open => Answered
>
>Luc Scholtès proposed the following answer:
>Hi,
>
>Yes, there is an easier way. Juste define the material parameters in
>such lines:
>
>def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson = ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
>def wallMat(): return JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))
>
>Second, make sure that you have these lines in your engines list:
>
>...
>InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
>InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
> [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
> [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
>),
>...
>
>with intR a value that you previously determined to ensure a
>coordination number K equal to the value you want for the packing you
>use (K=10 for sandstone).
>
>Then, the key parameter is: cohesiveTresholdIteration which defines the
>iteration up to which you want to create bonds between particles (the
>first one here).
>
>Thus, just by runnning
>
>O.step()
>
>the bonds will be created.
>
>After, you may want to reset the interaction radius to 1 by typing
>
>SSgeom.interactionDetectionFactor=-1.
>Saabb.aabbEnlargeFactor=-1.
>
>You can check that the coordination number K corresponds to the one you
>want by:
>
>numSSlinks=0
>numCohesivelinks=0
>for i in O.interactions:
> if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
> numSSlinks+=1
> if i.phys.isCohesive :
> numCohesivelinks+=1
>print "| K=", 2.0*numCohesivelinks/numSpheres
>
>and then run your simulation (triaxial loading).
>
>Luc
>
>--
>If this answers your question, please go to the following page to let us
>know that it is solved:
>https://answers.launchpad.net/yade/+question/556907/+confirm?answer_id=0
>
>If you still need help, you can reply to this email or go to the
>following page to enter your feedback:
>https://answers.launchpad.net/yade/+question/556907
>
>You received this question notification because you asked the question.

I guess that using the labels, you can get what you want. Here is a suggestion that actually I never used myself since I prepare my dense samples with a dedicated separate script (I then import these dense samples directly in my simulations, e.g. for triaxial deviatoric testing):

For the isotropic compaction phase, define default values in engines such as interactions take place between strictly contacting particles ( aabbEnlargeFactor=1, interactionDetectionFactor=1) and particles interact without cohesion (cohesiveTresholdIteration=0).

Then, when the confined state is reached (I guess you have a test for checking that), just use one step where you define the bonds at the right interaction range with the labels:

SSgeom.interactionDetectionFactor=intR
Saabb.aabbEnlargeFactor=intR
interactionPhys.cohesiveTresholdIteration=O.iter+1

O.step()

Then initialize the interaction range (if you want):

SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

Then run the deviatoric phase.

Let me know if it solves your problem.

Luc

>1. I can't visualize broken interactions in this way. At the end of the test, the interactions that showed >"isBroken" only make up about 5%

I noticed in JointedCohesiveFrictionalPM.cpp, that the "isBroken" flag is only tripped for bonds that fail according to the shear failure criteria (assuming neverErase=false). For interactions that fail according to the tensile failure criteria, the isBroken flag is only tripped if we use "neverErase =True".

This might explain why the isBroken != !isCohesive at the end of a triaxial test on a fully cohesive specimen.

Luc, maybe you can comment on why this is?

There actually is a sound justification to this behavior.

A shear bond failure is not a loss of interaction: the cohesive properties disappear but the interaction keeps existing (as long as there is geometrical overlap), we can still ask YADE for the interaction properties such as isBroken.

On the other hand, with a tensile failure and neverErase = False, the interaction immediately disappears: the tensile bond failure is equivalent to a loss of interaction. And how would you ask YADE the properties of something that does not exist ?

Jerome

Thanks, Jérôme, for the explanation.

Am I right to say there are no physical implications in the model when using neverErase=True? It is simply a flag that tells yade to keep track of old interactions?

I assume this will increase the storage of the model, does it reduce the speed of the simulation? Is yade recycling through these old interactions every time step?

I apologize if this is a thread derail, we can move this to a new thread if necessary.

I'm afraid it doesn't work quite well on my try. I wrote:

O.engines:

...

Bo1_Sphere_Aabb(label = 'bo1s'),...
Ig2_Sphere_Sphere_ScGeom(label = 'ig2ss'),...
Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=0, label = 'ip2jj'),...

...
triax,
newton,
PyRunner(..., label = 'checker'),

.......

def set_interaction():

bo1s.aabbEnlargeFactor = intRadius

ig2ss.interactionDetectionFactor = intRadius

ip2jj.cohesiveTresholdIteration = O.iter + 2

checker.iterPeriod = 1

checker.command = 'check_cohesion()'

def check_cohesion():

bo1s.aabbEnlargeFactor = -1

ig2ss.interactionDetectionFactor = -1

for i in O.interactions:

if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):

if i.phys.FnMax > 0:

num_coh += 1

And the num_coh showed 0.

--

At 2017-03-14 00:52:55, "Luc Scholtès" <email address hidden> wrote:
>Your question #556907 on Yade changed:
>https://answers.launchpad.net/yade/+question/556907
>
> Status: Open => Answered
>
>Luc Scholtès proposed the following answer:
>I guess that using the labels, you can get what you want. Here is a
>suggestion that actually I never used myself since I prepare my dense
>samples with a dedicated separate script (I then import these dense
>samples directly in my simulations, e.g. for triaxial deviatoric
>testing):
>
>For the isotropic compaction phase, define default values in engines
>such as interactions take place between strictly contacting particles (
>aabbEnlargeFactor=1, interactionDetectionFactor=1) and particles
>interact without cohesion (cohesiveTresholdIteration=0).
>
>Then, when the confined state is reached (I guess you have a test for
>checking that), just use one step where you define the bonds at the
>right interaction range with the labels:
>
>SSgeom.interactionDetectionFactor=intR
>Saabb.aabbEnlargeFactor=intR
>interactionPhys.cohesiveTresholdIteration=O.iter+1
>
>O.step()
>
>Then initialize the interaction range (if you want):
>
>SSgeom.interactionDetectionFactor=-1.
>Saabb.aabbEnlargeFactor=-1.
>
>Then run the deviatoric phase.
>
>Let me know if it solves your problem.
>
>Luc
>
>--
>If this answers your question, please go to the following page to let us
>know that it is solved:
>https://answers.launchpad.net/yade/+question/556907/+confirm?answer_id=2
>
>If you still need help, you can reply to this email or go to the
>following page to enter your feedback:
>https://answers.launchpad.net/yade/+question/556907
>
>You received this question notification because you asked the question.

Sorry Weijy but I got a bit lost...

It should not be so difficult to make this work.

Could you please send a working script so that I can have a look?

Luc

ps: Robert, the neverErase flag was introduced for using the DFNFlow engine and keep record of broken interactions when simulating hydraulic fracturing. Your question might need a dedicated thread.

>no physical implications in the model when using neverErase=True?

Yes correct. In practice it should be used combined with another model (e.g. capillary force model) which will effectively delete the interactions when needed. Else it is a memory leak as you suggest, but it may help bean counting.

There is a primitive implementation of callbacks in yade, they would be used for triggering the type of counter you need for particular events (creating/deleting interactions for instance) but I don't think they are functional at the moment.

Alternatively, you could implement the counting in the Law2 functor itself whenever something happen - an easy and quick solution.

Bruno

Thank you for the advice, Bruno.

Hi Luc, here is my script:

from yade import pack,export

compFricDegree = 30
confiningPressure = 1e4

O.materials.append(JCFpmMat(young=5e6, poisson=0.5, frictionAngle = radians(compFricDegree), cohesion = 1e3, tensileStrength = 1e2, density=2770, label='spheres'))
O.materials.append(FrictMat(young=5e6, poisson=0.5, frictionAngle = 0, density=0, label='walls'))
walls=aabbWalls([Vector3(0,0,0),Vector3(0.10,0.2,0.10)],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
sp.makeCloud(Vector3(0,0,0),Vector3(0.10,0.2,0.10),rRelFuzz = 0.2,num = 500,porosity = 0.95,seed = 1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
triax=TriaxialStressController(maxMultiplier=1.+5e-4,finalMaxMultiplier=1.+3e-5,thickness = 0,stressMask = 7,internalCompaction = True)
triax.goal1=triax.goal2=triax.goal3=-confiningPressure

O.engines=[
    ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(label = 'bo1s'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(label = 'ig2ss'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(label = 'ip2jj')],
  [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
 ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    NewtonIntegrator(),
    PyRunner(command = 'set_friction()', iterPeriod = 1000, label = 'checker'),
]

def set_friction():
    global compFricDegree
    print '*************************************************'
    print 'current step is ',O.iter
    print 'unbalanced force:',unbalancedForce(),' mean stress: ',triax.meanStress
    if unbalancedForce()<0.01 and abs(-confiningPressure-triax.meanStress)/confiningPressure<0.001:
        checker.command = 'set_interaction()'
        checker.ietrPeriod = 1

def set_interaction():
    num_intr = 0
    bo1s.aabbEnlargeFactor = 1.2
    ig2ss.interactionDetectionFactor = 1.2
    for i in O.interactions:
        num_intr = num_intr + 1
    print '*************************************************'
    print 'current step is ',O.iter
    print 'total interactions before setting intR is ', num_intr
    ip2jj.cohesiveTresholdIteration = O.iter + 1
    checker.iterPeriod = 1
    checker.command = 'check_cohesion()'

def check_cohesion():
    num_intr = 0
    num_coh = 0
    bo1s.aabbEnlargeFactor = -1
    ig2ss.interactionDetectionFactor = -1
    for i in O.interactions:
        if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
            num_intr += 1
            if i.phys.FnMax > 0:
                num_coh +=1
    print '*************************************************'
    print 'current step is ',O.iter
    print 'total interactions after setting intR is ', num_intr
    print 'number of cohesive interactions is ',num_coh
    O.pause()

Hi,

I would be surprised this is a working script as requested by Luc...
See the typo "checker.ietrPeriod = 1" in set_friction() definition..

Note also that it's kind of sad to use loops (in the set_interaction() definition) to count interactions in Yade (or elements of any Python list, from a more general point of view).
We fortunately have built-in functions for such purpose:
https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.InteractionContainer.countReal

Thank s for reminding, Jerome.

But I still don't get where is the mistake. I wrote "checker.ietrPeriod = 1" in set_friction() definition because it should(as I think) set the enlargeFactor on the very step and set cohesive interactions till the next step. And "checker.ietrPeriod = 1" ensures next step it comes to check_cohesion() and show the result of setting intRadius and cohesion. As I understand it should have the same effect as O.step(). I must got something wrong...

ps: Even I set "checker.ietrPeriod = 2" in set_friction(), result still goes the same.

I was mostly concerned by "ietrPeriod" which does not exist as an attribute for such an engine, as opposed to "iterPeriod" [*]....

[*] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.PyRunner.iterPeriod

That is awkward.. I copied a wrong script somehow. This should be the right one.
And a simple question: InteractionContainer.countReal does not return a value. How should I use it?

from yade import pack,export

confiningPressure = 1e4
O.materials.append(JCFpmMat(young=5e6, poisson=0.5, frictionAngle = radians(30), cohesion = 1e3, tensileStrength = 1e2, density=2770, label='spheres'))
O.materials.append(FrictMat(young=5e6, poisson=0.5, frictionAngle = 0, density=0, label='walls'))
walls=aabbWalls([Vector3(0,0,0),Vector3(0.10,0.2,0.10)],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
sp.makeCloud(Vector3(0,0,0),Vector3(0.10,0.2,0.10),rRelFuzz = 0.2,num = 500,porosity = 0.95,seed = 1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
triax=TriaxialStressController(maxMultiplier=1.+5e-4,finalMaxMultiplier=1.+3e-5,thickness = 0,stressMask = 7,internalCompaction = True)
triax.goal1=triax.goal2=triax.goal3=-confiningPressure

O.engines=[
    ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(label = 'bo1s'),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(label = 'ig2ss'),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(label = 'ip2jj')],
  [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
 ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    NewtonIntegrator(),
    PyRunner(command = 'set_interaction()', iterPeriod = 1000, label = 'checker'),
]

def set_interaction():
    print '*************************************************'
    print 'current step is ',O.iter
    print 'unbalanced force:',unbalancedForce(),' mean stress: ',triax.meanStress
    if unbalancedForce()<0.01 and abs(-confiningPressure-triax.meanStress)/confiningPressure<0.001:
        num_intr = 0
        bo1s.aabbEnlargeFactor = 1.2
        ig2ss.interactionDetectionFactor = 1.2
        for i in O.interactions:
            num_intr = num_intr + 1
        print '*************************************************'
        print 'current step is ',O.iter
        print 'total interactions before setting intR is ', num_intr
        ip2jj.cohesiveTresholdIteration = O.iter + 1
        checker.iterPeriod = 1
        checker.command = 'check_cohesion()'

def check_cohesion():
    num_intr = 0
    num_coh = 0
    bo1s.aabbEnlargeFactor = -1
    ig2ss.interactionDetectionFactor = -1
    for i in O.interactions:
        if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
            num_intr += 1
            if i.phys.FnMax > 0:
                num_coh +=1
    print '*************************************************'
    print 'current step is ',O.iter
    print 'total interactions after setting intR is ', num_intr
    print 'number of cohesive interactions is ',num_coh
    O.pause()

Thanks Luc Scholtès, that solved my question.