Bonding particles with JCFpm yields unexpected forces

Asked by David Alber on 2020-11-09

Dear Yade community,
being new to yade I lack some basic understanding about the constitutitve law ‘Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM’ and related JCFpmMat, Ip2_JCFpmMat_JCFpmMat_JCFpmPhys classes:

Problem Formulation:
Question 635871 [1] states that the JCFpm law is an implementation of Potyondy & Cundalls ‘A bonded particle model for rock’ [3]. I would therefore expect a Hertz Mindlin contact force [2. eq: 33-36], updated by a force proportional to a set bond normal (and shear) stiffness and the cross section A [3. eq: 14], as proposed by Potyondy & Cundall.
The following minimal test example with two spheres, each customly set ‘onJoint=True’ and the interaction set ‘isOnJoint=True’, shows that the normal force is proportional to: 2*young*(r1*r2/(r1+r2))*relativeOverlap. The stiffness kn is not calculated with the JCFpmMat parameter jointNormStiffness as expected from the source code JointedCohesiveFrictionalPM.ccp - line: 560-580.

Minimal example:
           jointCohesion=1.e8,jointFrictionAngle=radians(10.), label='Friction')
#Creating Spheres
s1 = utils.sphere([0.,0.,0.],radius=1,color=(1,0,1),material=mat,fixed=True)
s2 = utils.sphere([0.,0.,2.5],radius=1,color=(1,0,1),material=mat)
s2.state.vel = (0,0,-1.)
O.dt = .05*utils.PWaveTimeStep()
bo1s.aabbEnlargeFactor = 1.
ig2ss.interactionDetectionFactor = 1.
#Set sphere state and interactions manually 'onJoint'
for i in O.interactions:
   O.bodies[i.id1].state.joint += 1
   O.bodies[i.id2].state.joint += 1
   O.bodies[i.id1].state.onJoint, O.bodies[i.id2].state.onJoint = True, True
   i.phys.isOnJoint = True
#Check if force is related to jointNormalStiffness:,True)

i = O.interactions[0,1]
print('isOnJoint = ',i.phys.isOnJoint)
print('normalStifness=',, ' normalForce=',i.phys.normalForce, 'ref force=', (i.geom.penetrationDepth-i.phys.initD)*mat.young)

>>> isOnJoint = ', True, normalStifness=, 30000000000.0, normalForce=, Vector3(0,0,6063592.219396341), ref force=, 6063592.219396341

1) Is it necessary to custom set every interaction ‘isOnJoint’ when forces according to jointNormalStiffness / jointShearStiffness are desired?
2) How does one achieve forces that are guided by the cross-section A and a defined jointNormalStiffness?
3) Similar to the Mindlin Physics Model, is it possible to integrate viscous damping with a coefficient of restitution en

[1]: []
[2]: []
[3]: []

Question information

English Edit question
Yade Edit question
No assignee Edit question
Solved by:
Luc Scholtès
Last query:
Last reply:

This question was reopened

Luc Scholtès (luc) said : #1

Hi David,

1) no Hertz Mindlin contact force in JCFPM: stiffnesses are constant (linear force displacement law). The way these stiffnesses are computed is different for isOnJoint and not isOnJoint interaction though:

- for not isOnJoint interactions: kn is computed as in CundallStrack as a harmonic average of the particle elastic moduli such as kn=2* E1*R1*E2*R2 / (E1*R1 + E2*R2) and ks=ks = 2*E1*R1*v1*E2*R2*v2 / (E1*R1*v1 + E2*R2*v2)

- for isOnJoint interactions: kn and ks are directly assigned the values define by the jointNormal/ShearStiffness arguments such as ki = (jki1 + jki2)/2*contactPhysics->crossSection

2) Is it necessary to custom set every interaction ‘isOnJoint’ when forces according to jointNormalStiffness / jointShearStiffness are desired?

-> not sure to understand your inquiry here... I think I answered in the previous comment but: what exactly do you want to achieve?

3) How does one achieve forces that are guided by the cross-section A and a defined jointNormalStiffness?

-> I may not understand your questions but I think it is the case (see previous answer). You could also define the contact stiffness based on the surface when the interactions are set (after the first timestep of your simulation). You would have to do that for both the material attribute (body.mat) and the interaction attribute (interaction.phys) through python coding (no need to modify the sources).

4) Similar to the Mindlin Physics Model, is it possible to integrate viscous damping with a coefficient of restitution en (Ip2_FrictMat_FrictMat_MindlinPhys(en=en))?

-> yes, that's possible but you'll need to modify the sources if you want to combine JCFPM features with Mindlin Features (viscous damping + joint like interactions).

I hope that's clear


David Alber (davidalber) said : #2

Dear Luc,

I would have expected that. However, even if manually setting the interaction 'isOnJoint', the force between the two spheres is still calculated from the stiffness of the young module. This is demonstrated in the output of the minimal example:
Output: 'isOnjoint' = True; 'kn' = 3e10 (which is the young module and NOT jointNormalStiffness*phys->crossSection).
Is that related to the joint normal vector?

We want to compress a packing of spheres. A local non-viscous damping shall be represented by parallel bonds between every interaction to contribute to the energy loss in the system [1]. Hence, we seek for joints with set jointNormalStiffness between every interaction. This shall be achieved by the often cited formalism of Potyondy & Cundall 'A bonded particle model for rock' [2]. From the YADE launchpad question 635871 [3] I figured that the JCFpm related classes are an implementation of this formalism and would therefore be the right constitutive law.
Is that right?

3), 4) OK

[1]: []
[2]: []
[3]: []

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


1) After a quick look, problem seems to be that your interaction has already be assigned its mechanical properties (like, in the first O.step, when you declared it as onJoint.

Computation of these mechanical properties happen by default in YADE only once: at interaction creation (they usually do not change during a simulation => no need to recompute things which have already been computed).

Then, your manual changes on that interaction do not have the intended effect / are incomplete. With such a workflow, you also have to directly assign yourself an adequate kn value.

(Classical YADE workflow for the use of joint stiffness bodies parameters goes through examples/jointedCohesiveFrictionalPM/ script)

Luc Scholtès (luc) said : #4


For 1), as suggested by Jerome, if you want to modify an existing interaction stiffness, you need to add

in your loop

#Set sphere state and interactions manually 'onJoint'

That's true for all properties of an existing interaction. The way you did it only changed the properties on the particles, which will be taken into consideration the next time an interaction is created with these particles.

For 2), I am not sure to see the point of using a bond as a damping tool... For instance, you could compact the assembly with a classic frictional law using a high value of the non viscous damping coefficient (defined in Newton integrator)? Now, assuming that you really want to add cohesion between particles to compact the assembly , why define jointed interactions to do that? You can simply define cohesive bonds without the attribute "onJoint"... The onJoint attribute only makes sense if you want to reorientate the contacts according one or several predefined discontinuity planes (cf. smooth contact logic proposed in [1,2]). If that's not the case, just define cohesion and/or tensile strength between particles without changing the contact plane orientation (without playing on the onJoint attribute).



David Alber (davidalber) said : #5

@Jerome Duriez:
Thanks for your answer. Having a closer look at the example of the workflow, provided the answers I needed. It seems that not only the body property 'onJoint' needs to be true, but also the joint normal vectors need to point in the right directions. If both conditions are fulfilled the stiffness kn is calculated from the jointNormalStiffness.

@Luc Scholtes:
>just define cohesion and/or tensile strength between particles without changing the contact plane orientation (without playing on the onJoint attribute)<
Indeed the aim is to 'glue' particles together with cohesive bonds and changing the contact plane orientation is not a necessity.
What constitutive law would you recommend to achieve these cohesive bonds?
Are you referring to the Law2_ScGeom6D_CohFrictPhys_CohesionMoment law with the corresponding CohFrictMat material class?
As far as I understand, the parameter normalCohesion/ shearCohesion define a maximum tensile strength for the bond, which breaks if the force exceeds this value. Right?

Best Luc Scholtès (luc) said : #6

Hi David,

"What constitutive law would you recommend to achieve these cohesive bonds?"

You can use the JCFPM. All interactions (onJoint or notOnJoint) can be defined with or without cohesion/tensile strength.

The tensile strength defines a maximum admissible force in tension (normal direction to the contact).
The cohesion defines a maximum admissible force in shear (tangential direction to the contact).

Please find below an example for simulating a uniaxial compression test on an "intact" sample (no joints) with JCFPM.



# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-
from yade import pack, plot


#### packing (previously constructed)

#### Simulation Control
rate=-0.01 #deformation rate
iterMax=10000 # maximum number of iterations
saveVTK=2000 # saving output files for paraview

#### Material microproperties
intR=1.1 # allows near neighbour interaction (can be adjusted for every packing)
DENS=2500 # could be adapted to match material density: dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) -> packing porosity as to be computed?

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

#### create the specimen
O.bodies.append(pack.regularHexa(pred,radius=0.025,gap=0.,material=sphereMat)) # up to you to define another sample here, e.g., with randomDensePack or anything else.

for o in O.bodies:
 if isinstance(o.shape,Sphere):
   if o.shape.radius>Rmax:

print 'nbSpheres=',nbSpheres,' | Rmean=',Rmean

#### boundary condition (see utils.uniaxialTestFeatures

################# ENGINES DEFINED HERE

 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5, defaultDt=utils.PWaveTimeStep()),

################# RECORDER DEFINED HERE

def recorder():

# if you want to plot during simulation

################# PREPROCESSING

#### manage interaction detection factor during the first timestep and then set default interaction range ((cf. A DEM model for soft and hard rock, Scholtes & Donze, JMPS 2013))
### initializes the interaction detection factor

#### coordination number verification
for i in O.interactions:
    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
    if i.phys.isCohesive :

print "nblinks=", numSSlinks, " | nbCohesivelinks=", numCohesivelinks, " || Kcohesive=", 2.0*numCohesivelinks/nbSpheres


David Alber (davidalber) said : #7

Dear Luc Scholtes and YADE community,
This helps a lot. Thanks for the great support!

David Alber (davidalber) said : #8

Thanks Luc Scholtès, that solved my question.