Fix the Coulomb's criterion in linear, Hertz, and Hertz-Mindlin with limit force in DEM

[acdaigneault]

Description

The Coulomb's criterion was wrong in the particle-particle contact force for Hertz-Mindlin with limit force, Hertz and Linear.
The normal force norm was explicitly positive when doing normal_force.norm(), making the Coulomb's criterion always positive even if the particles are in repulsion, so in slidling.

Solution

Apply the the same way the Hertz-Mindlin with limit overlap is calculated.

Testing

Tests with low restitution coefficient may give negative coulomb criterion, this is why some tests have changed.

Documentation

Miscellaneous (will be removed when merged)

Checklist (will be removed when merged)

See this page for more information about the pull request process.

Code related list:

• All in-code documentation related to this PR is up to date (Doxygen format)
• Lethe documentation is up to date
• Fix has unit test(s) (preferred) or application test(s), and restart files are in the generator folder
• The branch is rebased onto master
• Changelog (CHANGELOG.md) is up to date
• Code is indented with indent-all and .prm files (examples and tests) with prm-indent

Pull request related list:

• Labels are applied
• There are at least 2 reviewers (or 1 if small feature) excluding the responsible for the merge
• If this PR closes an issue or is related to a project, it is linked in the "Projects" or "Development" section
• If the fix is temporary, an issue is opened
• The PR description is cleaned and ready for merge

[acdaigneault]

Heerreeee

[OGaboriault]

I'm not sure about the std::fabs ...

[OGaboriault]

wait, normal_force_norm can be negative??

normal_spring_constant is positive.
normal_overlap is positive.

I don't think the std::fabs is needed

[acdaigneault]

The constant damping is negative and this case is possible
|normal_spring_constant * normal_overlap| < |normal_damping_constant * normal_relative_velocity_value|

[acdaigneault]

No?

[OGaboriault]

yes the damping can be negative, but I don't thing we had an error.

With your implementation, we coulg go from friction to sliding back to friction right when the particle have near zero overlap but are moving appart.

I don't think this is physical.

[acdaigneault]

What do you mean with my implementation or what part?
This is the way Hertz-Mindlin limit overlap was implemented, but I added the fabs for the calculation of the coulomb limit.
Like, neighter of the model is wrong, I don't see how PP forces Hertz Limit Overlap and Limit Force culd be different up to here.
About the friction/sliding/friction, I will have to think about it tomorrow Zzz

[blaisb]

Ok so there's a slight different between the two, but the good approach I believe is actually the following:

    const double normal_force_value =
      normal_spring_constant * normal_overlap +
      normal_damping_constant * normal_relative_velocity_value;
normal_force = normal_force_value* normal_unit_vector;
    double coulomb_threshold =
      this->effective_coefficient_of_friction[vec_particle_type_index(
        particle_one_type, particle_two_type)] *normal_force_norm;

Notice that I renamed normal_force_norm to normal_force_value. Norm is a bad choice of a name here because a norm implies that the value is strictly positive. I am to blame for that, it's my fault.
The inconsistancy pointed out by @acdaigneault is a good one @OGaboriault. Neither of you are wrong, let me explain.

When the coefficient of restitution is low (<0.95) the Hertz model allow for an attractive force at the end of a contact. We know this is not physical, but that's the reality of it. It happens when the relative velocity is very large.
Consequently, the normal_force_value can be negative. However, it's alright that it is negative. However, it should not be used as an absolute value when establishing the Coulomb limit. Indeed, when the force is purely attractive (thus when the value above is negative), there is no friction (it makes no sense to have friction when you are essentially in traction and not in compression between the surfaces).

Consequently, the implementation that is bugged is the other one. AKA we should not be calculating the normal force tensor and taking it's norm, because that leads to the issue pointed out by @OGaboriault . What we should be doing is calculating the normal_force_value and then using this value to establish the Coulomb threshold, but without using any absolute value.

The implementation that was bugged is actually the calculate_hertz_mindlin_limit_force_contact which uses the norm and not the limit overlap. The latter is the one we should always use anyway so one could question if we really need to keep the limit force anyway.

Hope this clears things out. So in theory this PR could be closed and left as-is or can be remodified to fix the bug in the other implementation (the calculate_hertz_mindlin_limit_force_contact )

[acdaigneault]

Ok! Nice that we discuss that all together, It makes more sense that way.
I'm working on doing something to avoid to have the same calculation a some places, especially to avoid to kind of mistake. It is not the fault of anyone and it is actually always a good thing to find and fix them :)

What I can do is to fix the other forces in that PR so the changes on the tests will be easier to track or I can continue with my refactoring PR.
What do you perfer?

[OGaboriault]

By fix, you are talking about the name of the variable?

[acdaigneault]

The norm is explicitly calculated, forcing the positive value.
Yes the name of the variable is also wrong, but I'm not taling about it

[acdaigneault]

Yes the name have to be changed in any case

[blaisb]

@acdaigneault can you add a CHANGELOG entry then I would merge

[OGaboriault]

coulomb_threshold could be set to const.

Otherwise this is good to go.