[BUG] Kink in stress and strain fields at the boundaries

[danielwolff1]

Bug Description: What's going wrong?

I'm running a simple torsion test case where I fix a rectangular block at one end and prescribe a Dirichlet displacement at the opposite end which results in a torsion by 150 degrees. I'll attach my input script at the end of this issue.

When post processing the solution in paraview and evaluating the nodal stresses / strains along straight lines over the cross-sectional area, I noticed that the solution exhibits unphysical kinks, which after comparing with the element stresses / strains always seems to coincide with the element boundary (see attached screenshots).

Diagonal line at $y=2$, Cauchy stress component $\sigma_{xy}$

Another idea of mine was that something weird happens because the line lies diagonally in the domain and starts and ends in corners. However, even if I change the location of the line, the problem persists:

Horizontal line at $y=2$, Cauchy stress component $\sigma_{xx}$

It's also not a problem of a single field, I can look at different stress components and see the same behavior:

Horizontal line at $y=2$, Cauchy stress component $\sigma_{yz}$

Originally the line was located at $y=2$, which is the center of my domain. Even if I move it more towards one of the edges (e.g., $y=3$) the problem still persists:

Horizontal line at $y=3$, Cauchy stress component $\sigma_{yz}$

The mesh in the example that I attach is very coarse, but I refined the mesh and the behavior was always the same.

Expectations

I would have expected stress / strain curves without a kink at the edges of the computational domain. Since I am not a 4C developer, any insights are appreciated, e.g., if you have observed similar behavior in your simulations already and give a me a hint what else I could try out.

Possible Solution

Steps to Reproduce

• Input file: 📎 block_torsion.yaml
• Mesh: 📎 block_torsion_rename_ext_to_exo.txt

Your Environment

I was running my simulations on the IMCS cluster charon, using SHA c6d4500 from October 30th but we first noticed this behavior earlier this year in summer.

Related Issues and Pull Requests

Interested Parties

[amgebauer]

It might has something to do how we compute the stresses at the nodes and at the elements. On elements, we compute the simple average of all stresses at the Gauss points resulting in a single stress value for each element. For nodal stresses, we extrapolate the stress tensor using the shape functions and then do the averaging at the nodes (with the adjacent elements). Both approaches are physically not consistent and one should be aware of the implications!

What you could do is to write the stresses at the Gauss points. Then you can extrapolate the quantity to the nodes onto a discontinuous mesh (maybe with a consistent method) so that there will not be any averaging. This would be very interesting since it would also show you whether you have large differences on the nodes between the adjacent elements which would indicate locking in your case.

Note, that with our approach of averaging at elements and nodes, we have a hard time actually detecting locking from the stress output (both nodal and element stresses).