How to Dynamically Switch Contact Stiffness in pair_style granular/hertz Based on Overlap Threshold in LAMMPS

LAMMPS LAMMPS General Discussion
,
21

Hi jtclemm

Many thanks for your reponse, i have tried to use two wall to compress two particles this time, and it works, I can see clearly two piecewise stiffness behaviour, many thanks for your help and suggestions.
it just as you suggested, that its easy to modify to hertz under the pair_style granular, but would you mind share the idea whether it is possible to also have the piecewise inter-particle frictional coefficient.
saying before the threshold, we have the E_a and u_a (frictional coefficient), after the threshold, we have the E_b and u_b.
should we try to do modification on tangential part?

many thanks for your time and help
hope you have a nice day

best regards

deyun

22

Glad those tips helped!

One certainly could modify the frictional component and I don’t think it would be a significant challenge to concoct some functional form.

However, whether or not that is a physically reasonable idea I couldn’t say. That’d depend a lot on the exact details of the system you want to simulate and the science you hope to study. Topics that are beyond the scope of the LAMMPS forum.

23

Dear jtclemnn

Many thanks for your reposne, i have also successfully modified this piecewise stiffness for Hertz pair. Many thanks for your suggestions.
for the frictional component, we want to simulate coated particles, where rigid sand particle was covered by soft coating. Their modulus and frictional coefficient are different due to the experimental observation.

for the modification of the frictional? would you mind share sone idea?
i am thinking whether I can directly modify this gran_sub_mod_tangential.cpp.
I am thinking whether I can add a “if else” to check the overlap between particles at the very beginining before calculating the tangential forces. like
"void GranSubModTangentialMindlin::coeffs_to_local()
{
k = coeffs[0];
xt = coeffs[1];
mu1 = coeffs[2];
mu2 = coeffs[3];
delta_switch = coeffs[4]

if (gm->delta >= delta_switch) {

Emod = Emod2;
poiss = poiss2;
mu = mu2;

} else {

Emod = Emod1;
poiss = poiss1;
mu = mu1;

}
"
would you mind share the idea whether this is reasonable idea or not?

24

No, this would not be a good strategy. The method coeffs_to_local() is called during the setup phase at the start of a run. gm->delta would not be defined.

The least complicated method might be to define a new frictional tangential model that swaps between elastic properties as needed when forces are calculated based on the current value of gm->delta.

25

Dear jtclemm,

Thank you again for your kind response and helpful suggestions. Your clarification about gm->delta not being available during the setup phase makes perfect sense.

I hope you don’t mind me asking a further question that I’ve been struggling to understand.

In the function GranSubModTangentialMindlin::coeffs_to_local(), I noticed the following block of code:

if (k == -1) {
  if (!gm->normal_model->get_material_properties())
    error->all(FLERR, "Must either specify tangential stiffness or material properties...");

  double Emod = gm->normal_model->get_emod();
  double poiss = gm->normal_model->get_poiss();
}

It seems that Emod and Poiss are obtained from the normal contact model during this setup stage. May I ask — does this mean that while gm->delta is not available at this point, get_emod() and get_poiss() can still be used to access material properties from the normal model?

The reason I ask is that I’ve implemented a customized hertz/piecewise normal contact model, where Emod1, Emod2, Poiss1, and Poiss2 are selected based on a threshold delta_switch, but I did not define emod or poiss in this model. Still, the following command runs without any error:

pair_coeff * * hertz/piecewise ... tangential mindlin_rescale/force NULL 1.0 1.0

This left me a bit puzzled, and I would really appreciate your insights on why does this not raise an error, even though emod and poiss are not defined in my normal model? besides, In such a case, are any default values used for Emod and Poiss in the tangential component, or are they simply ignored? moreover, If I hope to allow tangential stiffness and friction coefficient to vary based on overlap (gm->delta), would it be acceptable to implement all logic in calculate_forces() rather than coeffs_to_local()?

Apologies if these questions are too basic — I’m still learning the internal design of LAMMPS granular sub-models and am trying to follow best practices.

Thanks again for your time and guidance.

best regards

deyun

26

Not at all, happy to help clarify the code.

The normal models calculate/define Emod and poiss during its own coeffs_to_local() method and the submodels are processed in order of normal → damping → tangential → rolling/twisting → heatflow (if I recall correctly). So the normal model will always have these variables defined.

If a normal model defines Emod and poiss, it needs to set a flag material_properties in the constructor. E.g. the hertz/material model sets it but hooke does not. I’m guessing your command runs because this flag is incorrectly turned on so LAMMPS thinks the variables Emod and poiss are defined so the tangential model is therefore grabs whatever it finds in memory.

I would recommend you either:

  1. Turn that flag off in your new model and manually specify a stiffness for the tangential model
  2. Pick one of set of elastic properties and save them to the original Emod and poiss variables so the tangential model grabs valid values.

Which option is better depends on the physics of your model.

Yes, that is correct. The coeffs_to_local() method is simply used to read/process/save the parameters of your contact model. Then, all other logic/math for forces/torques should be implemented in calculate_forces().

27

Dear jtclemm

I sincerely appreciate your help and valuable suggestions. Thanks to your guidance, I have initially successfully modified the LAMMPS code to achieve my goals. I will then do some verification for the modified code. I couldn’t have done it without your insights and patience.

Wishing you a wonderful day ahead

best regards

28

Dear jtclemm,

My apologies for bothering you again. While verifying my code, I encountered a rather odd issue when trying to restart a simulation.

The modified code works fine when I first run the simulation on the HPC. I use the command mpirun -np 48 ./lmp_mpi -in try1k.in, and everything runs as expected. I used the modified pair_style:

pair_coeff * * hertz/piecewise 3.5e8 7e10 0.25 0.25 1.0 0.005 tangential mindlin_rescale/force 3.2e8 6.4e10 1.0 0.001 0.001 0.005

Here, I set the values as follows:

  • Emod1 = 3.5e8
  • Emod2 = 7e10
  • Poiss1 = 0.25
  • Poiss2 = 0.25
  • Damp = 1.0
  • Delta_switch = 0.005

For the tangential part, following your suggestions, I specified:

  • k1 = 3.2e8
  • k2 = 6.4e10
  • xt = 1.0
  • mu1 = 0.001
  • mu2 = 0.001
  • Delta_switch = 0.005

However, when I stop the simulation and try to restart it, I always encounter the error: “Illegal Mindlin tangential model”. Upon checking the source code, I found that the error occurs when any of the following conditions are met:

if (k < 0.0 || xt < 0.0 || mu < 0.0) 
    error->all(FLERR, "Illegal Mindlin tangential model");

This is puzzling because all the k, xt, and mu values I provided are definitely larger than 0.

Then I tested the code on my own PC using ./lmp_mpi -in try1k.in, and encountered an unusual situation. Initially, the restart worked fine, but after stopping and trying to restart again, I randomly get the “Illegal Mindlin tangential model” error. Sometimes the restart works, sometimes it doesn’t. It’s inconsistent, and I’m unsure what could be causing it.

Additionally, when I ran mpirun -np 48 ./lmp_mpi -in try10k.in, I encountered a new error:

ERROR: Domain too large for neighbor bins.

This issue makes things more complicated. I would appreciate any insights or suggestions you might have regarding what’s really happening here. I attached all related files for your reference.
Thank you in advance for your time and help!

Best regards,
gran_sub_mod_custom.cpp (1.4 KB)
gran_sub_mod_custom.h (849 Bytes)
gran_sub_mod_tangential.cpp (14.0 KB)
gran_sub_mod_tangential.h (4.9 KB)
restart_filetry1k.52000000 (8.0 MB)
Particle_coarse.lj (3.6 MB)

try1k.in (4.8 KB)

29

So looking at your edits, it looks like you actually define k1, k2, mu1, and mu2 in gran_sub_mod_tangential.cpp lines 271-276. So the original variables, k and mu, are undefined and LAMMPS is just accessing whatever values randomly sit in memory (an invalid read). This likely explains the inconsistency you see.

You would want to update the variables in the error checks to your new variables, assign values to the original variables (e.g. k = k1), or something equivalent.

This likely implies something is wrong with your choice of simulation parameters. Maybe your pair coefficients or your barostat: 11.2. Errors and warnings details — LAMMPS documentation? I recommend visualizing your system during the lead up to this problem and/or turning off the barostat temporarily to start narrowing down what physically is going wrong.

Hope this helps.