Normal and Tangential force in granular pair style

Naghmeh · August 14, 2024, 9:59pm

Hello,

We are using the granular pair-style with the following setup:

pair_style granular
pair_coeff * * hertz/material $E ${ncdc} ${nu} tangential mindlin_rescale/force NULL ${tcdr} $f rolling sds ${kroll} ${rrdc} ${froll} twisting ${twistmodel}

For capturing the final geometry and topology of the system, we have defined:

compute final_geometry particles property/local patom1 patom2 cutoff radius
compute final_topology particles pair/local force p4 p10 p11 p12 cutoff radius

The outputs are dumped using:

dump dump_topo particles local 1 ${fname_write_dump_topo} c_final_geometry[*] c_final_topology[*]

It is known that in the in granular pair style, force represents the magnitude of the normal force acting between the pair of atoms, and P4 denotes the magnitude of the tangential force. According to the tangential mindlin_rescale/force option, the magnitude of the tangential force should be limited to μFn, where μ is the tangential friction coefficient and Fn is the normal force magnitude. Based on the final_topology definition, final_topology[1] should correspond to Fn, and final_topology[2] to Ft. Using a LAMMPS version dated 23 June 2022 gives results with Ft<Fn. However using a LAMMPS version of 28 March 2023 gives Ft>Fn.

Is it possible that the order of reporting Ft and Fn in final_topology has been swapped in the LAMMPS version from 28 March 2023? If not, what could be the potential issue?

Thanks

jtclemm · August 15, 2024, 6:15pm

Thank you for reporting the issue. Pair granular was recently overhauled and, in doing so, the output for the normal force accidentally omitted the normal damping force. This was corrected earlier this year: Coefficient of restitution based damping in granular models by dhairyaiitb · Pull Request #4068 · lammps/lammps · GitHub.

Does this rectify the error you see?

FWIW, this is the line that was updated:

github.com

lammps/lammps/blob/aa0b6c47c20cca07a7b2c63479ed7dd998094d3b/src/GRANULAR/pair_granular.cpp#L772


      
          model->vi = v[i];
          model->vj = v[j];
          model->omegai = omega[i];
          model->omegaj = omega[j];
          model->history = history;
          
          model->calculate_forces();
          
          // apply forces & torques
          // Calculate normal component, normalized by r
          fforce = model->Fntot * model->rinv;
          
          // set single_extra quantities
          svector[0] = model->fs[0];
          svector[1] = model->fs[1];
          svector[2] = model->fs[2];
          svector[3] = MathExtra::len3(model->fs);
          svector[4] = model->fr[0];
          svector[5] = model->fr[1];
          svector[6] = model->fr[2];
          svector[7] = MathExtra::len3(model->fr);

Naghmeh · August 15, 2024, 7:05pm

Thank you so much for the information.
Does this imply that we need to re-run all the simulations using a newer (or perhaps older) version of LAMMPS? Since we’ve already invested a lot of time in running these long simulations, I’m looking for a way to avoid unnecessary repetition.

jtclemm · August 15, 2024, 9:24pm

Without knowing the specific simulations/existing data/analysis, I cannot comment on what actions you might take. With that being said, this bug does not affect trajectories and only affects the reported normal force from this one command.

Depending on what data you output and what you need, you could use post processing to calculate damping forces and add them as required. Moving to another version is also an option. But if you have restart files for which you want to maintain compatibility, you could also just update this one method in the file and recompile. Either way, you may want to first confirm this change corrects the reported issue in your system.