# Bonding particles with JCFpm yields unexpected forces

Asked by David Alber on 2020-11-09

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:
tensileStrength=1.e8,cohesion=1.e8,jointNormalStiffness=2.e8,jointShearStiffness=2.e8,
#Creating Spheres
s2.state.vel = (0,0,-1.)
O.bodies.append([s1,s2])
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.5,label='bo1s'),]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.5,label='ig2ss'),],
[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(label='interactionPhys')],
[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(label='interactionLaw')]
),
NewtonIntegrator(damping=0.),
]
O.dt = .05*utils.PWaveTimeStep()
O.step()
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:
O.run(10,True)

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

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

Questions:
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
(Ip2_FrictMat_FrictMat_MindlinPhys(en=en))?

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Luc Scholtès
Solved:
2020-11-18
Last query:
2020-11-18
2020-11-18

## This question was reopened

 Luc Scholtès (luc) said on 2020-11-10: #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

Luc

 David Alber (davidalber) said on 2020-11-12: #2

Dear Luc,

1)
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?

2)
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

 Jérôme Duriez (jduriez) said on 2020-11-12: #3

Hi

1) After a quick look, problem seems to be that your interaction has already be assigned its mechanical properties (like phys.kn), 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/identifBis.py script)

 Luc Scholtès (luc) said on 2020-11-12: #4

Alright,

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

i.phys.kn=valueYouWant

#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).

Luc

 David Alber (davidalber) said on 2020-11-17: #5

@Jerome Duriez:
1)
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:
2)
>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?

 Luc Scholtès (luc) said on 2020-11-18: #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.

Luc

---

# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-

################# SIMULATIONS DEFINED HERE

#### packing (previously constructed)
OUT='compressionTest_JCFPM'

#### 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?
YOUNG=20e9
FRICT=7
ALPHA=0.1
TENS=1e6
COH=1e6

#### material definition

#### create the specimen
pred=pack.inCylinder((0,0,0),(0,1,0),0.25)
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.

R=0
Rmax=0
nbSpheres=0.
for o in O.bodies:
if isinstance(o.shape,Sphere):
nbSpheres+=1
Rmean=R/nbSpheres

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

#### boundary condition (see utils.uniaxialTestFeatures
bb=utils.uniaxialTestFeatures()
negIds,posIds,longerAxis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

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

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(recordCracks=True,Key=OUT,label='interactionLaw')]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5, defaultDt=utils.PWaveTimeStep()),
NewtonIntegrator(damping=0.4,label='newton'),
PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
]

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

def recorder():
'eps':strainer.strain,
'sigma':strainer.avgStress,
'tc':interactionLaw.nbTensCracks,
'sc':interactionLaw.nbShearCracks,
'te':interactionLaw.totalTensCracksE,
'se':interactionLaw.totalShearCracksE,
'unbF':utils.unbalancedForce()})
plot.saveDataTxt(OUT)

# if you want to plot during simulation
plot.plots={'i':('sigma')}
plot.plot()

################# 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))
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

#### 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 :

################# SIMULATION REALLY STARTS HERE
O.run(iterMax)

 David Alber (davidalber) said on 2020-11-18: #7

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

 David Alber (davidalber) said on 2020-11-18: #8

Thanks Luc Scholtès, that solved my question.