Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fixed stress visualization for elastic materials #4375

Merged
merged 7 commits into from
Nov 30, 2023
Merged

fixed stress visualization for elastic materials #4375

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
e5302bd
fixed stress visualization for elastic materials
bobmyhill Oct 1, 2021
ad35c1e
corrected sign convention, added test
bobmyhill Nov 10, 2023
3ea11c6
added changelog
bobmyhill Nov 10, 2023
3391ba8
added shear stress postprocessor
bobmyhill Nov 10, 2023
0cfc1ab
added stress output
bobmyhill Nov 10, 2023
3e94a8e
added viscoelastic stress visualization tests
bobmyhill Nov 11, 2023
a7bf368
allow for non-zero initialized values
bobmyhill Nov 29, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions doc/modules/changes/20231110_bobmyhill
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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>();
bobmyhill marked this conversation as resolved.
Show resolved Hide resolved

// 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
Open in desktop
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
Loading