Hi all,
We (@mibillen) have been looking at the principal stress visualization outputs in our 2D visco-elastic-plastic subduction models and were surprised by the principal stress orientations calculated in the slab. We ran a 2D visco-plastic subduction model for comparison, and there was significant (90 degrees in some areas) differences in principal stress orientations between a visco-elastic-model and a viscous model (with the viscous model showing the physically more expected stress orientations in the slab). The visco-elastic orientations seemed wrong, and we’re surprised, but if you’re not comparing viscous to visco-elastic models you might not notice. I have now run some small test models to illustrate what we were seeing and what seems like it might be incorrect.
I modified the bending beam box model to suspend the beam as a dense block of rock from the top of box. We would expect this beam to be in down-dip tension as gravity is pulling it downward. I ran two versions of this model that are identical except for the “Enable elasticity” formulation being set to true in one model, and false in the other. ( VEtopbeam.prm (6.6 KB), Vtopbeam.prm (6.6 KB))
These outputs show the results of running the two versions of the beam model. You can see on the left, in the viscoelastic model, that the sigma_1 direction indicates the beam is in down-dip compression, contrary to what we physically would expect. On the right, in the viscous model, the sigma_1 direction indicates the beam is in down-dip tension, as we would physically expect.
In aspect/source/postproccess/visualization/principal_stress.cc, we saw the note on compressive stress being positive followed by the flag for incorporating elasticity:
I tested changing all of the “+=” signs within this elastic part to “-=”, recompiled, and then re-ran the above models. The new “principal_stress_direction_1” orientations are rotated and now show down-dip tension in the hanging beam both in the viscoelastic and in the viscous model, as expected:
@jbnaliboff @anne-glerum We think this is the correct orientation and that this sign in incorporating the elasticity was just incorrect, unless there is something physical I am misunderstanding between the viscous and viscoelastic behavior?
Additionally, this is just in the post processing, but I was not sure where to find within the solvers how the stress tensor was being handled sign wise while solving the stokes system?
