Kelvin-Voigt Boundary Implementation for #692

[cgeudeker]

Describe the PR
Kelvin-Voigt boundary implementation for #692

Related Issues/PRs
This PR aims to implement the Kelvin-Voigt boundary, or more specifically the spring-dashpot system, as a boundary conditions. This utilizes linear elastic material properties, as added by PR #705, as well as json file input parameters (shown in the Additional context section) to calculate traction values for the following equations:

Spring coefficients kp and ks are calculated according to the following equations:

Other relevant to this PR, are the parameters delta, a and b. Each of these will be user input values under the boundary conditions as shown in the additional context section below. a and b are optional input parameters, which if not specified will be equivalent to a value of 1. However, the following relation must be satisfied:

where h_min is taken to be the cell height.

Additional context
input.json file is currently expected to look like the following

"boundary_conditions": {
...
   "absorbing_boundary": [
              {
                  "nset_id" : #,
                  "dir" : #,
                  "delta": #,
                  "a": #,
                  "b": #
              }
          ],
...
},

The above implementation is based on the work of this paper.

[kks32]

Add unit tests and please fix as per suggestions

[kks32]

When do we map this information?

[cgeudeker]

I wasn't able to get the unit test done today, but I wanted to make sure that some of the other edits I was able to make today were what was expected since I was slightly unsure, particularly with regards to apply_absorbing_boundary. I also still wasn't quite sure where to map the new material information pwave_velocity and swave_velocity onto the nodal properties. I thought I had understood it to be in mpm_explicit.tcc but I couldn't find anything similar.

[codecov]

Codecov Report

Merging #710 (4cbc60a) into develop (adfe6f1) will decrease coverage by 0.12%.
The diff coverage is 86.44%.

@@             Coverage Diff             @@
##           develop     #710      +/-   ##
===========================================
- Coverage    96.75%   96.63%   -0.12%     
===========================================
  Files          130      131       +1     
  Lines        25932    26248     +316     
===========================================
+ Hits         25090    25363     +273     
- Misses         842      885      +43     

Impacted Files	Coverage Δ
include/data_types.h	66.67% <ø> (ø)
include/loads_bcs/constraints.h	100.00% <ø> (ø)
include/mesh.h	100.00% <ø> (ø)
include/node.h	100.00% <ø> (ø)
include/node_base.h	100.00% <ø> (ø)
include/particles/particle.h	100.00% <ø> (ø)
include/particles/particle_base.h	100.00% <ø> (ø)
include/solvers/mpm_base.h	50.00% <ø> (ø)
include/solvers/mpm_scheme/mpm_scheme.tcc	60.92% <0.00%> (-7.03%)	⬇️
include/solvers/mpm_base.tcc	73.63% <23.81%> (-3.41%)	⬇️
... and 9 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 adfe6f1...4cbc60a. Read the comment docs.

[jgiven100]

@cgeudeker Would it work to follow along with nodal_velocity_constraints() in mpm_base.tcc here?

Within mpm_base it calls

constraints_->assign_nodal_velocity_constraints()

then in constraints.tcc it calls

mesh_->nodes(set_id)->assign_velocity_constraint()

at which point the velocity constaints are saved as nodal properties. For the absorbing boundary case, we would need to define 2 new variables in node.h to save the pwave and swave values.

[cgeudeker]

A few notes about the last commit.

I see that I need to put together a test case that will go through the absorbing constraint file, so I'll look at what was done for the others and follow along one of those. I'll probably commit this later today.

There is an issue that resurfaced while finalizing the test case for the mapping of wave velocities to nodes. When 2 particles of the same material were mapped to the nodes, the wave velocities doubled. I thought this would be fixed by having the update_property boolean argument set to false, but the same result occurred. So if there is any advice about how to best solve this it would be much appreciated.

Finally, some of the sections are poorly commented, so that will be fixed in the next couple of commits, but I wanted to make sure that the function was working as desired.

[jgiven100]

@kks32 2 mains issues have come up. Guidance on these would be appreciated.

1: Nodal location affects how to go from traction -> force

Problem

The thesis provides an equation of traction for the absorbing boundary. We need to apply this as a nodal force. Based on the nodes location (corner, side, or bottom), the different components of the traction must be multiplied by either half or full cell lengths to become force. Below are the 3 cases to consider.

Proposed solution

Add "position" : to input JSON to indicate if the specific node set is for the model's "corner", "side", or "bottom". For example:

      "absorbing_constraints" : [
        { 
          "nset_id" : 0,
          "dir" : 1,
          "delta" : 0.25,
          "h_min" : 0.001,
          "a" : 1.0,
          "b" : 1.0,
          "position" : "corner"
        }
      ]

Then update apply_absorbing_constraint(...) to handle the 3 possible cases.

2: How to deal with material points displacing beyond the absorbing boundary

Problem

If delta is very small, then wave energy at the boundary will be reflected instead of absorbed. So delta must be somewhat large to be effective. But this also means that material points will likely displace beyond the absorbing boundary. If the mesh also extends beyond the absorbing boundary, the model will continue to run, but this breaks conservation of mass within the sub-domain defined by the absorbing boundary.

1. Is extending the mesh beyond the absorbing boundary an accetable work-around?
2. Is conservation of mass within the subdomain defined by the absorbing boundary necessary?
3. If extending mesh beyond the absorbing boundary is ok, would that affect the area of influence for going from traction to force (problem above)? (i.e., displacement and velocity will be mapped to the node form inside and outside the boundary, but the traction would still be only applied to the partial cell within the boundary.)

side note: I don't think Q3 is important for geotechnical applications, but rather for keeping the code general. Absorbing boundaries are likely to only be applied to the sides of models with lateral damping being a priority. Area of influence for lateral damping will be along vertical cell edges which wouldn't be affected by material points outside the absorbing boundary.

EDITS:

• (Nov. 23, 2021): cross out wrong scheme

[jgiven100]

@kks32 I've been able to locally implement the proposed solution described above (adding "position" to boundary inputs). I'll wait to hear from you if this is an acceptable solution before pushing to the feature branch.

[kks32]

@jgiven100 That's a good solution. Thanks!!

[jgiven100]

@kks Seems like there might have been a Clang update? Clang is failing here. Seems to have been passing before.

[kks32]

@kks Seems like there might have been a Clang update? Clang is failing here. Seems to have been passing before.

I thought it was cppcheck not clang:

include/node.tcc:404:37: warning: Found suspicious operator ',' [constStatement]
            absorbing_force_ << dir0, dir1;
                                    ^
include/node.tcc:409:37: warning: Found suspicious operator ',' [constStatement]
            absorbing_force_ << dir0, dir1;

You can turn-off clang-format temporarily using

    // clang-format off
    ...code...
    // clang-format on

[jgiven100]

@kks32 @cgeudeker

Feature Updates

• I've locally updated node.tcc to handle 3D models.

Verification

• With a 3D column, I've been able to reproduce 1D p-wave results (from Al-Kafaji thesis).
• I'll work on upward moving shear wave next. Once reproducing the 2D results, I'll push changes here (next few days).

Next Steps

• Current 2D benchmark Kelvin-Voigt Boundary Benchmark Case mpm-benchmarks#20 needs some slight updates due to changes in input file naming convention.
• Add identical 3D benchmark to the above open bencmark PR.
• Open mpm-doc PR.

Edits

(24 Nov. 2021) - make next steps into check boxes

[jgiven100]

@kks @cgeudeker
This is ready for a re-review.

I'm happy with unit testing coverage for 18 of 21 changed files. Expections are mpm_base.tcc, mpm_scheme.tcc, and mpm_explicit.tcc. Looking at similar existing functions for the exception files, I don't see too much existing testing... This leads me to believe that the areas with low coverage in this update is ok.