reference for cohesion model in LIGGGHTS?

Submitted by ericparteli on Wed, 08/14/2013 - 16:29

Dear Christoph, dear Andreas,
Could you please provide some reference(s) that use(s) or derive(s) the equations of the cohesion models used in LIGGGHTS (skrj and skrj2 equations; see below)?
Thanks so much!
A simple linear cohesion model can be activated by setting cohesion = 'sjkr' or 'sjkr2' (simplified Johnson-Kendall-Roberts model). If two particle are in contact, it adds an additional normal force tending to maintain the contact, which writes

F = k A,

where A is the particle contact area and k is the cohesion energy density in J/m3. For sjkr, the sphere-sphere contact area is calculated as (

A = Pi/4 * ((dist-Ri-Rj)*(dist+Ri-Rj)*(dist-Ri+Rj)*(dist+Ri+Rj) )/(dist*dist)

where dist is the distance between the particle centers. For sjkr2, the sphere-sphere contact area is calculated as

A = 2*Pi * delta_n * (2R*)

If you are using the linear cohesion model, you must also define the cohesion energy density:

fix id all property/global cohesionEnergyDensity peratomtypepair n_atomtypes value_11 value_12 .. value_21 value_22 .. .
(value_ij=value for the cohesion energy density (in Energy/Length3 units) between atom type i and j; n_atomtypes is the number of atom types you want to use in your simulation)

Maryam | Wed, 09/11/2013 - 04:22

Hi Eric,

Thanks for the very good question! As I had difficulty rapping my head around this simplified cohesion model and estimating the cohesion energy density (without tuning it to experimental data), I referred to this article which gives a review of available cohesion models, including JKR model:

But I am no less confused now; as far as I understood from the derivations in the above paper, the cohesion force is proportional to a^1.5 while in the LIGGGHTS formulation it is a function of a^2 (a being the contact radius), so I can’t estimate the cohesion energy density from other parameters in the original JKR formulation. I might go ahead and modify the formulation in the code.

I’d appreciate any clarifying input.


ericparteli | Tue, 09/17/2013 - 09:45

Hi, Maryam,

Thanks for your comment.

Well what I did was to implement in LIGGGHTS a cohesion model which I called sjkr3, in which I take for the cohesion force the original expression derived from the JKR model: F = 4 sqrt( PI * gamma * a^3 * Y_eff ), with a = sqrt (R_eff * delta_n) standing for the contact radius.

You can find this expression for example in the second term of Eq. (5) of this paper:
X. L. Deng and R. N. Davé, Dynamic simulation of particle packing influenced by size, aspect ratio and surface energy, Granular matter 15, 401-415 (2013).

Note that the quantity gamma there is the surface energy in J/m^2.

However, it would be very helpful if someone could explain to us where the expressions of the current public release of LIGGGHTS (sjkr and sjkr2) come from!


Maryam | Thu, 09/19/2013 - 00:20

yes, the equation that I mentioned above also has the same form. Can you comment on the results that you got? Did you get to match them to any experimental data?


ericparteli | Thu, 09/26/2013 - 10:20

Hi, Maryam,
>> yes, the equation that I mentioned above also has the same form.
OK. I should note, in correction to my last post, that the contact radius should be calculated in consistence with the JKR result, see e.g. Eq.(35) of your reference, Barthel (2008). In fact the contact radius sqrt(R_eff * delta_n) is obtained from the cohesionless Hertz model...
>> Can you comment on the results that you got? Did you get to match them to any experimental data?
I am right now working in this direction, I hope to post good news soon!

rahulsoni | Fri, 05/22/2015 - 17:34

Hello Eric

The standard models, as mentioned in references, use surface energy (J/m2). However, LIGGGHTS uses an uncommon term energy density (J/m3), whose values are difficult to get. Can anyone modify the models as per the standard models shown in mentioned references.


Rahul Kumar Soni
Scientist, CSIR - IMMT, India

ckloss's picture

ckloss | Mon, 06/15/2015 - 14:58

Hi Rahul,

you will have to calibrate anyway - no matter if the input values are common or uncommon. You can also modify the cohesion model yourself - it's relatively easy


golovnya | Tue, 02/16/2016 - 11:44

Dear Christoph, dear Andreas,

The question is opened in 2013 and it still retaines topicality. Could you, please, explain does the application of sjkr and sjkr2 give correct results?
And please, provide some references for k (J/m3) or relation between k and gamma (J/m2), which is usually used.

Thank you in advance.

hrvig | Fri, 02/26/2016 - 11:24


For those who don't want to implement the JKR model in LIGGGHTS, I suggest that you compare the output of the two models. I compared the the sjkr model (F=k_sjkr*A) to the jkr model described in J. S. Marshall "Discrete-element modeling of particulate aerosol flows", Journal of Computational Physics 228, 1541-1561; 2009.

That way, you can specify a Hamaker constant and get the surface energy density (gamma) in J/m^2. Using that one you can fit the cohesion energy density (J/m^3) to get the same force as function of overlap.


golovnya | Mon, 06/27/2016 - 16:10

I am sorry for the late reply. Thank you very much for the reference and details. It is still up-to-date.

MiRa | Thu, 12/22/2016 - 12:05


I also found it hard to reasonably connect the implemented JKR cohesion models with the equations in the paper of Johnson et al. [1]. According to Equation 18 of this paper, the "apparent Hertz load" is a function of the inital load, which compresses both bodies.

I just did a quick and dirty research about this topic, but would it not be possible for the developers to implement the equation(s) from this publication? What was the reason for not using these original equations?

Best regards,

[1] K. L. Johnson, K. Kendall and A. D. Roberts (1971). Surface Energy and the Contact of Elastic Solids. Proceedings of the Royal Society of London. Series A, Mathematical and PhysicalSciences, Vol. 324, No. 1558, pp. 301-313

ckloss's picture

ckloss | Mon, 01/02/2017 - 13:26

Hi MiRa,

sure, any alternative approach to the SJKR model in LIGGGHTS should be pretty straight forward to implement. Feel free to contribute :-)

Best wishes

DavidGaube1 | Fri, 01/27/2017 - 21:15


I have a question about the sjkr2 cohesion model. The contact area is calculated as A = 2*Pi * delta_n * (2R*).
I can't understand where the second "2" in the equation comes from. When I derive that equation, I only get one "2".