Contact mechanics softening factor

[IvarStefansson]

In the contact mechanics equations, we currently compare tractions to c_num times (increments of) displacement jumps. C_num is by default computed as
c_num = shear modulus / (solid.residual_aperture() + solid.fracture_gap()) * SOFTENING FACTOR.

From a physical intuition, this resembles a Hooke's law, thus ensuring that the tractions and scaled displacements have similar magnitudes. However, the justification of the softening factor is truly heuristic ("making the fracture softer than the matrix " for the default value of 0.1). We should consider to either purge it or replace the shear modulus with the fracture's elastic modulus, i.e. the fracture stiffness parameter (which typically takes values roughly one order of magnitude smaller than the Lame parameters).

[IvarStefansson]

This is partly related to #957

[IvarStefansson]

I think we can solve this based on the now scaled version as well as some improved understanding of the problem. Let's discuss in person, @keileg and @jwboth.

[IvarStefansson]

Suggested action:

• Remove softening factor
• Write down mechanical equations with all scalings
• Improve documentation

[jwboth]

Some preliminary observations/conclusions - still WIP and not ready for discussion but as a reminder for later discussion:

• characteristic_contact_traction among the material parameters is not used in the code and can be removed.
• characteristic_displacement in exact sense does not enter longer the effective definition of c_num, as indicated in the issue description. When rescaling the contact equations with characteristic_contact_traction as implemented in the core, c_num has the interpretation c_num = shear modulus / L * SOFTENING FACTOR, where L is the max length of the geometry. Numerically, the performance of Newton's method remains depening on the choice of characteristic_displacement. Yet, not too strongly. One guess why, is conditioning of the contact condition. Also the characteristic displacement controls the angle of the disncontinuity in the derivative in the transformed coordinates.
• If characteristic_displacement is chosen on the order of the normal jump, it indeed is balancing the arguments of the normal contact condition, and choosing softeningfactor equals to 1, is thus natural. There is no guarantee that the tangential increment is of the same order though, which has to be kept in mind. Guessing the right sizes, may be possible in academic cases, but cannot be done though global rescaling. Thus, it is tempting to make a more thorough parameter study which varies the above parameters and investigate the sensitivity.

[IvarStefansson]

To your first point: It was intended for the case where you know a traction driving force and not the displacement (e.g. you have only zero Dirichlet values and nonzero Neumann). In that case, you might want to override the methods for characteristic displacement and traction, fetching the traction material parameter for the latter and estimating the former from traction based on Hooke's law (inverse of what is by default done in the traction method). Does this make sense?

[jwboth]

Thanks for the comment. Without a NOTE in the code, this is not obvious (In fact in previous run scripts, I have set values for both characteristic quantities although I thought I checked their use, but overlooked the lacking effect of changing the traction value) - it was not for us. We will add a note then.

[jwboth]

Alternative could be to add the rescaling you describe and make a if-else choice. If the characteristic traction != 1, one would invoke the added rescaling, otherwise, the default rescaling is used. If both values are uneqal one, a ValueError is thrown.

[mariusnevland]

Alternative could be to add the rescaling you describe and make a if-else choice. If the characteristic traction != 1, one would invoke the added rescaling, otherwise, the default rescaling is used. If both values are uneqal one, a ValueError is thrown.

Related to this point, would it be a good idea to include in the if-else condition that if both characteristic_displacement and characteristic_traction are equal to one, then these are the characteristic values that will be used? In this way, the default option will be to use the contact equations in their original, non-scaled form.

[mariusnevland]

Or at least I think there should be an obvious way to use the non-scaled contact equations if one wishes to do so. At the moment it is necessary to override the characteristic_contact_traction() method to achieve this, if I have understood correctly.

[jwboth]

Since the work on the rescaling did not only add rescaling, but also changed the meaning of cnum, I do not think it would be advisable to re-establish the old cnum when both parameters are chosen equal to 1. This would invoke a discontinuous behavior. On the other hand, my description (or having the two different scalings depending on the two different inputs) is also not continuous... We can discuss next time - I guess it depends on relevance and use cases. So far the philosopy was the user is responsible for adding the mixins he/she needs. The current indication of what mixin should be developed and added breaks with this philsophy. This is why I was suggesting to actually add the intended use.

[jwboth]

Or at least I think there should be an obvious way to use the non-scaled contact equations if one wishes to do so. At the moment it is necessary to override the characteristic_contact_traction() method to achieve this, if I have understood correctly.

One would have to be careful and overwrite both characteristic_contact_traction and characteristic_displacement.