Skip to content

Commit

Permalink
Merge pull request #4375 from bobmyhill/fix_vep_stress_viz
Browse files Browse the repository at this point in the history 
fixed stress visualization for elastic materials
  • Loading branch information
@gassmoeller
gassmoeller authored Nov 30, 2023
2 parents 49ff2a9 + a7bf368 commit e6eadff
Show file tree
Hide file tree
Showing 12 changed files with 48,525 additions and 45,047 deletions.
4 changes: 4 additions & 0 deletions doc/modules/changes/20231110_bobmyhill
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
Fixed: The stress and shear stress visualization postprocessors
now output the correct stresses when elasticity is enabled.
<br>
(Bob Myhill, Rebecca Fildes, 2023/11/10)
49 changes: 29 additions & 20 deletions source/postprocess/visualization/shear_stress.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,33 +67,42 @@ namespace aspect
// ...and use it to compute the stresses
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q];
const SymmetricTensor<2,dim> deviatoric_strain_rate
= (this->get_material_model().is_compressible()
?
strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor<dim>()
:
strain_rate);

const double eta = out.viscosities[q];

// Compressive stress is positive in geoscience applications
SymmetricTensor<2,dim> shear_stress = -2.*eta*deviatoric_strain_rate;

// Add elastic stresses if existent
// Compressive stress is negative by the sign convention
// used by the engineering community, and as input and used
// internally by ASPECT.
// Here, we change the sign of the stress to match the
// sign convention used by the geoscience community.
SymmetricTensor<2,dim> shear_stress;

// If elasticity is enabled, the deviatoric stress is stored
// in compositional fields, otherwise the deviatoric stress
// can be obtained from the viscosity and strain rate.
if (this->get_parameters().enable_elasticity == true)
{
shear_stress[0][0] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];
shear_stress[1][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];
shear_stress[0][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];
shear_stress[0][0] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];
shear_stress[1][1] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];
shear_stress[0][1] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];

if (dim == 3)
{
shear_stress[2][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];
shear_stress[0][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];
shear_stress[1][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];
shear_stress[2][2] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];
shear_stress[0][2] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];
shear_stress[1][2] = -in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];
}
}
else
{
const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q];
const SymmetricTensor<2,dim> deviatoric_strain_rate
= (this->get_material_model().is_compressible()
?
strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor<dim>()
:
strain_rate);

const double eta = out.viscosities[q];
shear_stress = -2. * eta * deviatoric_strain_rate;
}

for (unsigned int d=0; d<dim; ++d)
for (unsigned int e=0; e<dim; ++e)
Expand Down
52 changes: 30 additions & 22 deletions source/postprocess/visualization/stress.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,34 +67,42 @@ namespace aspect
// ...and use it to compute the stresses
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q];
const SymmetricTensor<2,dim> deviatoric_strain_rate
= (this->get_material_model().is_compressible()
?
strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor<dim>()
:
strain_rate);

const double eta = out.viscosities[q];

// Compressive stress is positive in geoscience applications
SymmetricTensor<2,dim> stress = -2.*eta*deviatoric_strain_rate +
in.pressure[q] * unit_symmetric_tensor<dim>();

// Add elastic stresses if existent
if (this->get_parameters().enable_elasticity == true)
// Compressive stress is negative by the sign convention
// used by the engineering community, and as input and used
// internally by ASPECT.
// Here, we change the sign of the stress to match the
// sign convention used by the geoscience community.
SymmetricTensor<2,dim> stress = in.pressure[q] * unit_symmetric_tensor<dim>();

// If elasticity is enabled, the deviatoric stress is stored
// in compositional fields, otherwise the deviatoric stress
// can be obtained from the viscosity and strain rate.
if (this->get_parameters().enable_elasticity)
{
stress[0][0] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];
stress[1][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];
stress[0][1] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];
stress[0][0] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xx")];
stress[1][1] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yy")];
stress[0][1] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xy")];

if (dim == 3)
{
stress[2][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];
stress[0][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];
stress[1][2] += in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];
stress[2][2] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_zz")];
stress[0][2] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_xz")];
stress[1][2] -= in.composition[q][this->introspection().compositional_index_for_name("ve_stress_yz")];
}
}
else
{
const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q];
const SymmetricTensor<2,dim> deviatoric_strain_rate
= (this->get_material_model().is_compressible()
?
strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor<dim>()
:
strain_rate);

const double eta = out.viscosities[q];
stress -= 2. * eta * deviatoric_strain_rate;
}

for (unsigned int d=0; d<dim; ++d)
for (unsigned int e=0; e<dim; ++e)
Expand Down
7 changes: 4 additions & 3 deletions tests/viscoelastic_stress_build-up.prm
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
#
# The analytical solution for viscoelastic stress build-up in
# an incompressible medium with a constant strain-rate is:
# simga_ij = 2 * edot_ii * eta * (1 - e^(-mu*t/eta)),
# sigma_ij = 2 * edot_ii * eta * (1 - e^(-mu*t/eta)),
# where sigma_ij is the elastic deviatoric stress tensor, edot_ij is
# the deviatoric strain-rate, eta is the background viscosity, mu is the
# shear modulus and t is time.
Expand All @@ -35,8 +35,8 @@
# in the compositional fields tracking viscoelastic stresses.
#
# For the current test, we prescribe initial values to the compositional fields
# that track the viscoelastic stresses that are equal to the eventual
# analytical solution for the viscoelastic stress.
# that track the viscoelastic stresses that are equal to the
# analytical solution for the viscoelastic stress after 250 Kyr.

include $ASPECT_SOURCE_DIR/benchmarks/viscoelastic_stress_build-up/viscoelastic_stress_build-up.prm

Expand Down Expand Up @@ -77,6 +77,7 @@ end
subsection Postprocess
subsection Visualization
set Time between graphical output = 1e3
set List of output variables = stress, shear stress, stress second invariant
set Output format = gnuplot
end
end
Loading

0 comments on commit e6eadff

Please sign in to comment.