Bonding particles with JCFpm yields unexpected forces

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

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

[1]: [https://www.sciencedirect.com/science/article/pii/S0921883118302061]
[2]: [https://www.sciencedirect.com/science/article/pii/S1365160904002874]
[3]: [https://answers.launchpad.net/yade/+question/635871]

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)

Alright,

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

i.phys.kn=valueYouWant

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

Luc

[1] https://www.sciencedirect.com/science/article/pii/S1365160910002169
[2] https://www.sciencedirect.com/science/article/pii/S1365160912000391?via%3Dihub

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

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

Thanks Luc Scholtès, that solved my question.