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

Jump to bottom

Fix elastic stress used in principal_stress postprocessor #5480

Merged
merged 1 commit into from
Nov 19, 2023
Merged

Fix elastic stress used in principal_stress postprocessor #5480

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
6 changes: 6 additions & 0 deletions doc/modules/changes/20231109_gassmoeller
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
Fixed: The 'principal stress' postprocessor used the wrong
sign for stresses in models with elastic deformation, leading
to principal stress directions that were rotated by 90 degrees.
This is fixed now.
<br>
(Rene Gassmoeller, Rebecca Fildes, 2023/11/09)
57 changes: 34 additions & 23 deletions source/postprocess/visualization/principal_stress.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -115,36 +115,39 @@ namespace aspect

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>();
SymmetricTensor<2,dim> stress;

// Add elastic stresses if existent
if (this->get_parameters().enable_elasticity == true)
if (this->get_parameters().enable_elasticity == false)
{
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")];
const SymmetricTensor<2,dim> deviatoric_strain_rate
= (this->get_material_model().is_compressible()
?
in.strain_rate[q] - 1./3. * trace(in.strain_rate[q]) * unit_symmetric_tensor<dim>()
:
in.strain_rate[q]);

// Compressive stress is positive in geoscience applications
stress = -2. * out.viscosities[q] * deviatoric_strain_rate;
}
else
{
// Compressive stress is positive in geoscience applications
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")];
}
}

if (use_deviatoric_stress == true)
stress -= 1./dim * trace(stress) * unit_symmetric_tensor<dim>();

if (use_deviatoric_stress == false)
stress += in.pressure[q] * unit_symmetric_tensor<dim>();

const std::array<std::pair<double, Tensor<1,dim>>, dim> principal_stresses_and_directions = eigenvectors(stress);

Expand Down Expand Up @@ -219,9 +222,17 @@ namespace aspect
ASPECT_REGISTER_VISUALIZATION_POSTPROCESSOR(PrincipalStress,
"principal stress",
"A visualization output object that outputs the "
"principal stress values and directions, i.e., the "
"principal stresses and directions, i.e., the "
"eigenvalues and eigenvectors of the stress tensor. "
"The postprocessor can either operate on the full "
"Wikipedia defines principal stresses as follows: "
"At every point in a stressed body there are at least "
"three planes, called principal planes, with normal "
"vectors, called principal directions, where the "
"corresponding stress vector is perpendicular to the "
"plane, and where there are no normal "
"shear stresses. The three stresses normal to these "
"principal planes are called principal stresses. "
"This postprocessor can either operate on the full "
"stress tensor or only on the deviatoric stress tensor, "
"depending on what run-time parameters are set."
"\n\n"
Expand Down
155 changes: 155 additions & 0 deletions tests/visco_elastic_top_beam.prm
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,155 @@
# Like the viscous_top_beam.prm test, but including elasticity.


# Global parameters
set Dimension = 2
set End time = 2e3
set Use years in output instead of seconds = true
set Output directory = output_ve

set CFL number = 0.5
set Maximum time step = 1e3

# Solver settings
subsection Solver parameters
subsection Stokes solver parameters
set Number of cheap Stokes solver steps = 2000
end
end

# Model geometry (7.5x5 km, 0.1 km spacing)
subsection Geometry model
set Model name = box

subsection Box
set X repetitions = 15
set Y repetitions = 10
set X extent = 7.5e3
set Y extent = 5e3
end
end

# Mesh refinement specifications
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 0
set Time steps between mesh refinement = 0
end

# Element types
subsection Discretization
set Temperature polynomial degree = 1
set Use locally conservative discretization = false
set Use discontinuous temperature discretization = false
set Use discontinuous composition discretization = true

subsection Stabilization parameters
set Use limiter for discontinuous composition solution = true
set Global composition maximum = 1.e11, 1.e11, 1.e11, 1.0
set Global composition minimum = -1.e11, -1.e11, -1.e11, 0.0
end
end

# Formulation classification
subsection Formulation
set Enable elasticity = true
end

# Velocity boundary conditions
subsection Boundary velocity model
set Zero velocity boundary indicators = top
set Tangential velocity boundary indicators = bottom, left, right
end

# Number and name of compositional fields
subsection Compositional fields
set Number of fields = 4
set Names of fields = ve_stress_xx, ve_stress_yy, ve_stress_xy, beam
end

# Spatial domain of different compositional fields
subsection Initial composition model
set Model name = function

subsection Function
set Variable names = x,y
set Function constants =
set Function expression = 0; 0; 0; if (x>=3.5e3 && x<=4.0e3 && y>=3.0e3, 1, 0)
end
end

# Composition boundary conditions
subsection Boundary composition model
set Fixed composition boundary indicators = bottom, top, right, left
set List of model names = initial composition
end

# Temperature boundary conditions
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top, left, right
set List of model names = box

subsection Box
set Bottom temperature = 293
set Left temperature = 293
set Right temperature = 293
set Top temperature = 293
end
end

# Temperature initial conditions
subsection Initial temperature model
set Model name = function

subsection Function
set Function expression = 293
end
end

# Material model
subsection Material model
set Model name = visco plastic

subsection Visco Plastic
set Densities = 2800, 2800, 2800, 2800, 6000
set Viscous flow law = dislocation
set Prefactors for dislocation creep = 5.e-19, 5.e-19, 5.e-19, 5.e-19, 5.e-25
set Stress exponents for dislocation creep = 1.0
set Activation energies for dislocation creep = 0.
set Activation volumes for dislocation creep = 0.
set Elastic shear moduli = 1.e11, 1.e11, 1.e11, 1.e11, 1.e10
set Fixed elastic time step = 2e3
set Use fixed elastic time step = true
set Viscosity averaging scheme = maximum composition
set Cohesions = 1.e50
end
end

# Gravity model
subsection Gravity model
set Model name = vertical

subsection Vertical
set Magnitude = 10.
end
end

# Post processing
subsection Postprocess
set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization

subsection Visualization
set List of output variables = material properties, strain rate, stress second invariant, principal stress
set Output format = gnuplot
set Time steps between graphical output = 1
set Interpolate output = true

subsection Principal stress
set Use deviatoric stress = true
end

subsection Material properties
set List of material properties = density, viscosity
end
end
end
50 changes: 50 additions & 0 deletions tests/visco_elastic_top_beam/screen-output
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@

Number of active cells: 150 (on 1 levels)
Number of degrees of freedom: 7,054 (1,302+176+176+1,350+1,350+1,350+1,350)

*** Timestep 0: t=0 years, dt=0 years
Solving temperature system... 0 iterations.
Skipping ve_stress_xx composition solve because RHS is zero.
Skipping ve_stress_yy composition solve because RHS is zero.
Skipping ve_stress_xy composition solve because RHS is zero.
Solving beam system ... 0 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 35+0 iterations.

Postprocessing:
RMS, max velocity: 0.000658 m/year, 0.002 m/year
Temperature min/avg/max: 293 K, 293 K, 293 K
Writing graphical output: output-visco_elastic_top_beam/solution/solution-00000

*** Timestep 1: t=1000 years, dt=1000 years
Solving temperature system... 0 iterations.
Solving ve_stress_xx system ... 2 iterations.
Solving ve_stress_yy system ... 2 iterations.
Solving ve_stress_xy system ... 3 iterations.
Solving beam system ... 2 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 27+0 iterations.

Postprocessing:
RMS, max velocity: 0.000543 m/year, 0.00167 m/year
Temperature min/avg/max: 293 K, 293 K, 293 K
Writing graphical output: output-visco_elastic_top_beam/solution/solution-00001

*** Timestep 2: t=2000 years, dt=1000 years
Solving temperature system... 0 iterations.
Solving ve_stress_xx system ... 2 iterations.
Solving ve_stress_yy system ... 2 iterations.
Solving ve_stress_xy system ... 2 iterations.
Solving beam system ... 2 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 26+0 iterations.

Postprocessing:
RMS, max velocity: 0.0005 m/year, 0.00148 m/year
Temperature min/avg/max: 293 K, 293 K, 293 K
Writing graphical output: output-visco_elastic_top_beam/solution/solution-00002

Termination requested by criterion: end time