Feature/contact

[thiagordonho]

Describe the PR Implementation of the contact-related functions in the Node and Particle classes and definition of the ContactFriction class. These entities were incorporated to MPMExplicit and are currently able to deal with frictional contact between two distinct bodies.

Related Issues/PRs This PR is related to RFC #636 .

Additional context This PR closely follows the design detail presented in RFC #636 . In fact, it comprehends the implementations of the contact-related functions to the Node and Particle classes. The functions in the particles are responsible for mapping the properties from the material points to the nodes (started with map_multimaterial) and updating the particle's kinematics using the acceleration and velocity from nodes (compute_contact_updated_position). The functions in the nodes are responsible for computing the additional nodal properties (compute_multimaterial_velocity, compute_multimaterial_normal_unit_vector, compute_multimaterial_relative_velocity, compute_contact_acceleration_velocity) that are required for resolving contact at contact nodes using apply_contact_mechanics.

To validate the contact algorithm implemented, two 2D benchmark problems were implemented: a block sliding down an inclined plane, and a cylinder rolling down an inclined plane. Both problems were modeled in the horizontal position with an inclined gravity applied to simulate the inclination of the plane. Both models also considered the Maximum Velocity Gradient, or MVG, as the normal computation type (see RFC #636 ). They were tested for different friction coefficients and angles of inclination. The results for the position of the center of mass were compared to that of the analytic solutions for rigid block and cylinder cases. Figures 1 and 2 show a schematic of the model from both problems and Figures 3 to 6 show the obtained results. As shown in Figure 2, the cylinder was modeled at 1 cell above the plane because the approximation condition would not handle the rolling properly if the cylinder is directly touching the plane. Other approximation conditions must be implemented to solve this issue.

Figure 1: Sliding block model in its initial position.

Figure 2: Rolling cylinder model in its initial position.

Figure 3: Horizontal position of the center of mass of the sliding block for a 30-degree inclination.

Figure 4: Horizontal position of the center of mass of the sliding block for a 45-degree inclination.

Figure 5: Horizontal position of the center of mass of the rolling cylinder for a 30-degree inclination.

Figure 6: Horizontal position of the center of mass of the rolling cylinder for a 45-degree inclination.

[codecov[bot]]

Codecov Report

Merging #716 (d0de0d7) into develop (86ba10e) will decrease coverage by 0.16%. The diff coverage is 95.56%.

@@             Coverage Diff             @@
##           develop     #716      +/-   ##
===========================================
- Coverage    96.78%   96.62%   -0.16%     
===========================================
  Files          130      131       +1     
  Lines        25899    26548     +649     
===========================================
+ Hits         25065    25650     +585     
- Misses         834      898      +64     

Impacted Files	Coverage Δ
include/mesh.h	100.00% <ø> (ø)
include/node_base.h	100.00% <ø> (ø)
include/particles/particle.h	100.00% <ø> (ø)
include/particles/particle_base.h	100.00% <ø> (ø)
include/contacts/contact.h	71.43% <66.67%> (-28.57%)	:arrow_down:
include/solvers/mpm_explicit.tcc	90.00% <70.00%> (-5.59%)	:arrow_down:
include/solvers/mpm_scheme/mpm_scheme.tcc	67.90% <83.33%> (-0.05%)	:arrow_down:
include/particles/particle.tcc	93.43% <87.36%> (-1.79%)	:arrow_down:
include/mesh.tcc	83.18% <90.91%> (+0.19%)	:arrow_up:
include/contacts/contact_friction.tcc	93.62% <91.67%> (-6.38%)	:arrow_down:
... and 15 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 86ba10e...d0de0d7. Read the comment docs.

[jgiven100]

@thiagordonho Thanks for all of the figures on this. For the cylinder cases, the limiting friction angle and coefficient are given as:

\phi_L = atan( sin( \alpha / 3 ) )
\mu_L = tan( atan( sin( \alpha / 3 ) ) )
\mu_L = sin( \alpha / 3 )

For the examples you show, this implies the following (please let me know if I'm thinking about this correctly):

45 Deg. Incline

\mu_L = sin( 45 / 3 )
\mu_L = 0.26

so \mu=0.2 is "sliding" while \mu=0.4 is "rolling".

Question

For the cylinder + 45 deg. incline, I think its good that the \mu = 0.2 case (red dash, Figure 6) aligns with the analytical solution like the sliding block outputs since in this case it is just sliding, like the blocks. I expect most contact scenarios we will model with CB-Geo will be a sliding mechanism, so great to catpure this so well!

But I'm a bit suprised the "rolling" cases are over-predicting displacement. I would have expected rolling to under-predict. Do you have any ideas why numerical outputs always appear to over-predict here?