← Back to team overview

yade-users team mailing list archive

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

 

New question #693898 on Yade:
https://answers.launchpad.net/yade/+question/693898

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:
mat=JCFpmMat(young=3.e10,poisson=0.25,density=5000.,frictionAngle=radians(10.),
           tensileStrength=1.e8,cohesion=1.e8,jointNormalStiffness=2.e8,jointShearStiffness=2.e8,
           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.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)

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

[1]: [https://answers.launchpad.net/yade/+question/635871]
[2]: [https://onlinelibrary.wiley.com/doi/epdf/10.1002/ceat.201200116]
[3]: [https://www.sciencedirect.com/science/article/pii/S1365160904002874]


-- 
You received this question notification because your team yade-users is
an answer contact for Yade.