Skip to content

Commit

Permalink
Merge pull request geodynamics#5480 from gassmoeller/fix_principal_st…
Browse files Browse the repository at this point in the history 
…ress_elastic

Fix elastic stress used in principal_stress postprocessor
  • Loading branch information
@bangerth
bangerth authored Nov 19, 2023
2 parents 0c68b8d + b8442ef commit 299a645
Show file tree
Hide file tree
Showing 13 changed files with 9,682 additions and 1,373 deletions.
6 changes: 6 additions & 0 deletions doc/modules/changes/20231109_gassmoeller
View file
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
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
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
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



Loading

0 comments on commit 299a645

Please sign in to comment.