How to add a new contact model

Bias's picture
Submitted by Bias on Wed, 07/04/2018 - 08:32

Hi,
I want to add a new contact model but I don't understand the file structure for liggghts. What I want to do is to add a full JKR model.
I was thinking then to make a copy of "normal_model_hertz.h", rename it to normal_model_jkr and work in this one.
What do I need to change to be able to use this model with a input file as "pair_style gran model jkr"?

I guess I should change the start of the file:
#ifdef NORMAL_MODEL
NORMAL_MODEL(HERTZ,hertz,3)
#else
#ifndef NORMAL_MODEL_HERTZ_H_
#define NORMAL_MODEL_HERTZ_H_

But what number should I put? I tried looking into "normal_model_base.h" but I feel a little rusty in C++ so I don't understand exactly. Are there any other files needed to be modified inorder to add a model?

What would be the difference if I added a jkr model as cohesion for example cohesion_model_jkr.h and use that with Hertz like
"pair_style gran model hertz cohesion jkr"?
Both normal_model_hertz.h and cohesion_model_jkr.h would then have the function surfacesIntersect but would the cohesion code overwrite the normal_model or how does that work?

Best regards
Tobias Eidevåg

Forums:

medvedeg | Wed, 07/04/2018 - 10:26

Hi Tobias,

you can use any unique integer number, not shared by any other normal force model.
If you implement a cohesion model inside a normal force model then the cohesion force is only active when particles do overlap. In some cohesion models, the cohesion force can be active even if the particles do not overlap. But this is not the case for sjkr cohesion model. I think you can safely implement your hybrid normal/cohesion model.

Alexander Podlozhnyuk

Bias's picture

Bias | Thu, 07/05/2018 - 08:03

Hi Alexander,
Thanks for the info. So for adding a new contact model its only needed to add a new file and that file will be detected automatically when compiling?

Another question: How is it best to save values for debugging a model? I would like to be able to add temporary parameters that are saved during a simulation for example with the "dump" command. Is that possible to do?
Or is there another way to debug models in a good way? I would for example as a test run 1 single sphere against a wall and then save the overlap, contact area and the contact force for every time step.

Best regards
Tobias Eidevåg

Bias's picture

Bias | Mon, 07/09/2018 - 18:14

Hi again,
I have another question related to this and to what Alexander wrote above.

So I want to add the JKR model and I did so by copying the normal_model_hertz.h to a new file normal_model_jkr.h. Now I can compile and run that model so that is great so far.
I see however that I don't get any negative overlap as I want to have. So I guess liggghts stop using the function surfacesIntersect when there is no overlap. I would however like it to run a little longer just until a certain threshold overlap.

How do I do that?

arnom's picture

arnom | Thu, 07/19/2018 - 09:27

Then the surfacesClose function will be used. Note that you need to call the function neighbor->register_contact_dist_factor(...) in your connectToProperties function in order to assure that non-touching contacts are handled well.

DCS team member & LIGGGHTS(R) core developer

Bias's picture

Bias | Wed, 08/15/2018 - 16:48

Hi,
I managed to implement the contact model as descired above so thanks all for the help on that.

One of the challenges with the JKR model I have now is that is very slow and I don't think I can run with many particles. I read an article by J.S Marshall "Discrete-element modeling of particulate aerosol flows" where they used a look up table for the simulation. So a normalized contact force and radius is precomputed and stored in a look up table as a function of displacement.

Can this be done in LIGGGHTS somehow?
I'm not exactly sure how to proceed and how to best store this look up table. Have something similar been done before? I was thinking if I could store a matrix variable that is only defined in the beginning and then interpolate values from that one based on the current displacement but I dont know where to define this variable so its only executed once for the simulation.

paul | Sat, 08/18/2018 - 20:58

> where to define this variable so its only executed once for the simulation.
Add a class member to your cohesionModel, and initialize its values in the constructor is my first guess.
But proceed w/ caution, you only have access to things that are defined before the respective pair_style / fix wall/gran* calls, as this is when this constructor is called.

If everything is known beforehand, consider using static constexpr at compile-time, although I don't know how much this will inflate the lmp_* binary size. I suspect this can be an issue due to the heavy use of templating, then a different class (containing the table) to be included could be the way to go.

khalifa | Wed, 05/01/2019 - 13:49

Dear Tobias,

I am trying to implement the full JKR model and I am having some doubts. Do you mind sharing your code/modifications?

I am specifically interested in how you solved the 4th power equation in order to compute the radius of contact from the overlap? In addition, did you manage to use a lookup table?

Your answers will be much appreciated. Thanks in advance.

Ali