# yade-users team mailing list archive

## Re: [Question #693898]: Bonding particles with JCFpm yields unexpected forces

```Question #693898 on Yade changed:

Description changed to:
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))?