Skip to content

Commit a5bc338

Browse files
committed
Add shear heating plugin
1 parent ca2225f commit a5bc338

File tree

3 files changed

+177
-113
lines changed

3 files changed

+177
-113
lines changed

cookbooks/anisotropic_viscosity/av_material.cc

Lines changed: 0 additions & 113 deletions
Original file line numberDiff line numberDiff line change
@@ -48,8 +48,6 @@
4848
#include <deal.II/base/geometry_info.h>
4949

5050
#include <aspect/material_model/simple.h>
51-
#include <aspect/heating_model/interface.h>
52-
#include <aspect/heating_model/shear_heating.h>
5351
#include <aspect/gravity_model/interface.h>
5452
#include <aspect/simulator/assemblers/stokes.h>
5553
#include <aspect/simulator_signals.h>
@@ -60,103 +58,6 @@
6058

6159
namespace aspect
6260
{
63-
namespace HeatingModel
64-
{
65-
template <int dim>
66-
class ShearHeatingAnisotropicViscosity : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
67-
{
68-
public:
69-
/**
70-
* Compute the heating model outputs for this class.
71-
*/
72-
void
73-
evaluate (const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
74-
const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
75-
HeatingModel::HeatingModelOutputs &heating_model_outputs) const override;
76-
77-
/**
78-
* Allow the heating model to attach additional material model outputs.
79-
*/
80-
void
81-
create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &material_model_outputs) const override;
82-
};
83-
84-
85-
86-
template <int dim>
87-
void
88-
ShearHeatingAnisotropicViscosity<dim>::
89-
evaluate (const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
90-
const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
91-
HeatingModel::HeatingModelOutputs &heating_model_outputs) const
92-
{
93-
Assert(heating_model_outputs.heating_source_terms.size() == material_model_inputs.position.size(),
94-
ExcMessage ("Heating outputs need to have the same number of entries as the material model inputs."));
95-
96-
// Some material models provide dislocation viscosities and boundary area work fractions
97-
// as additional material outputs. If they are attached, use them.
98-
const std::shared_ptr<const ShearHeatingOutputs<dim>> shear_heating_out
99-
= material_model_outputs.template get_additional_output_object<ShearHeatingOutputs<dim>>();
100-
101-
const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity
102-
= material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>();
103-
104-
for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q)
105-
{
106-
// If there is an anisotropic viscosity, use it to compute the correct stress
107-
const SymmetricTensor<2,dim> &directed_strain_rate = ((anisotropic_viscosity != nullptr)
108-
?
109-
anisotropic_viscosity->stress_strain_directors[q]
110-
* material_model_inputs.strain_rate[q]
111-
:
112-
material_model_inputs.strain_rate[q]);
113-
114-
const SymmetricTensor<2,dim> stress =
115-
2 * material_model_outputs.viscosities[q] *
116-
(this->get_material_model().is_compressible()
117-
?
118-
directed_strain_rate - 1./3. * trace(directed_strain_rate) * unit_symmetric_tensor<dim>()
119-
:
120-
directed_strain_rate);
121-
122-
const SymmetricTensor<2,dim> deviatoric_strain_rate =
123-
(this->get_material_model().is_compressible()
124-
?
125-
material_model_inputs.strain_rate[q]
126-
- 1./3. * trace(material_model_inputs.strain_rate[q]) * unit_symmetric_tensor<dim>()
127-
:
128-
material_model_inputs.strain_rate[q]);
129-
130-
heating_model_outputs.heating_source_terms[q] = stress * deviatoric_strain_rate;
131-
132-
// If shear heating work fractions are provided, reduce the
133-
// overall heating by this amount (which is assumed to be converted into other forms of energy)
134-
if (shear_heating_out != nullptr)
135-
heating_model_outputs.heating_source_terms[q] *= shear_heating_out->shear_heating_work_fractions[q];
136-
137-
heating_model_outputs.lhs_latent_heat_terms[q] = 0.0;
138-
}
139-
}
140-
141-
142-
143-
template <int dim>
144-
void
145-
ShearHeatingAnisotropicViscosity<dim>::
146-
create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &material_model_outputs) const
147-
{
148-
const unsigned int n_points = material_model_outputs.viscosities.size();
149-
150-
if (material_model_outputs.template has_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>() == false)
151-
{
152-
material_model_outputs.additional_outputs.push_back(
153-
std::make_unique<MaterialModel::AnisotropicViscosity<dim>> (n_points));
154-
}
155-
156-
this->get_material_model().create_additional_named_outputs(material_model_outputs);
157-
}
158-
}
159-
16061
namespace MaterialModel
16162
{
16263
// The AV material model calculates an anisotropic viscosity tensor from director vectors and the normal and shear
@@ -481,20 +382,6 @@ namespace aspect
481382
// explicit instantiations
482383
namespace aspect
483384
{
484-
namespace HeatingModel
485-
{
486-
ASPECT_REGISTER_HEATING_MODEL(ShearHeatingAnisotropicViscosity,
487-
"anisotropic shear heating",
488-
"Implementation of a standard model for shear heating. "
489-
"Adds the term: "
490-
"$ 2 \\eta \\left( \\varepsilon - \\frac{1}{3} \\text{tr} "
491-
"\\varepsilon \\mathbf 1 \\right) : \\left( \\varepsilon - \\frac{1}{3} "
492-
"\\text{tr} \\varepsilon \\mathbf 1 \\right)$ to the "
493-
"right-hand side of the temperature equation.")
494-
}
495-
496-
497-
498385
namespace MaterialModel
499386
{
500387
ASPECT_REGISTER_MATERIAL_MODEL(AV,
Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
/*
2+
Copyright (C) 2025 by the authors of the ASPECT code.
3+
4+
This file is part of ASPECT.
5+
6+
ASPECT is free software; you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation; either version 2, or (at your option)
9+
any later version.
10+
11+
ASPECT is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with ASPECT; see the file LICENSE. If not see
18+
<http://www.gnu.org/licenses/>.
19+
*/
20+
21+
#ifndef _aspect_heating_model_shear_heating_anisotropic_viscosity_h
22+
#define _aspect_heating_model_shear_heating_anisotropic_viscosity_h
23+
24+
#include <aspect/heating_model/interface.h>
25+
26+
#include <aspect/simulator_access.h>
27+
28+
namespace aspect
29+
{
30+
namespace HeatingModel
31+
{
32+
template <int dim>
33+
class ShearHeatingAnisotropicViscosity : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
34+
{
35+
public:
36+
/**
37+
* Compute the heating model outputs for this class.
38+
*/
39+
void
40+
evaluate (const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
41+
const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
42+
HeatingModel::HeatingModelOutputs &heating_model_outputs) const override;
43+
44+
/**
45+
* Allow the heating model to attach additional material model outputs.
46+
*/
47+
void
48+
create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &material_model_outputs) const override;
49+
};
50+
}
51+
}
52+
53+
#endif
Lines changed: 124 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,124 @@
1+
/*
2+
Copyright (C) 2015 - 2023 by the authors of the ASPECT code.
3+
4+
This file is part of ASPECT.
5+
6+
ASPECT is free software; you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation; either version 2, or (at your option)
9+
any later version.
10+
11+
ASPECT is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with ASPECT; see the file LICENSE. If not see
18+
<http://www.gnu.org/licenses/>.
19+
*/
20+
21+
#include <aspect/heating_model/shear_heating_anisotropic_viscosity.h>
22+
23+
#include <aspect/material_model/interface.h>
24+
#include <aspect/simulator/assemblers/stokes_anisotropic_viscosity.h>
25+
26+
#include <deal.II/base/symmetric_tensor.h>
27+
#include <deal.II/base/signaling_nan.h>
28+
29+
namespace aspect
30+
{
31+
namespace HeatingModel
32+
{
33+
template <int dim>
34+
void
35+
ShearHeatingAnisotropicViscosity<dim>::
36+
evaluate (const MaterialModel::MaterialModelInputs<dim> &material_model_inputs,
37+
const MaterialModel::MaterialModelOutputs<dim> &material_model_outputs,
38+
HeatingModel::HeatingModelOutputs &heating_model_outputs) const
39+
{
40+
Assert(heating_model_outputs.heating_source_terms.size() == material_model_inputs.position.size(),
41+
ExcMessage ("Heating outputs need to have the same number of entries as the material model inputs."));
42+
43+
// Some material models provide dislocation viscosities and boundary area work fractions
44+
// as additional material outputs. If they are attached, use them.
45+
const std::shared_ptr<const ShearHeatingOutputs<dim>> shear_heating_out
46+
= material_model_outputs.template get_additional_output_object<ShearHeatingOutputs<dim>>();
47+
48+
const std::shared_ptr<const MaterialModel::AnisotropicViscosity<dim>> anisotropic_viscosity
49+
= material_model_outputs.template get_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>();
50+
51+
for (unsigned int q=0; q<heating_model_outputs.heating_source_terms.size(); ++q)
52+
{
53+
// If there is an anisotropic viscosity, use it to compute the correct stress
54+
const SymmetricTensor<2,dim> &directed_strain_rate = ((anisotropic_viscosity != nullptr)
55+
?
56+
anisotropic_viscosity->stress_strain_directors[q]
57+
* material_model_inputs.strain_rate[q]
58+
:
59+
material_model_inputs.strain_rate[q]);
60+
61+
const SymmetricTensor<2,dim> stress =
62+
2 * material_model_outputs.viscosities[q] *
63+
(this->get_material_model().is_compressible()
64+
?
65+
directed_strain_rate - 1./3. * trace(directed_strain_rate) * unit_symmetric_tensor<dim>()
66+
:
67+
directed_strain_rate);
68+
69+
const SymmetricTensor<2,dim> deviatoric_strain_rate =
70+
(this->get_material_model().is_compressible()
71+
?
72+
material_model_inputs.strain_rate[q]
73+
- 1./3. * trace(material_model_inputs.strain_rate[q]) * unit_symmetric_tensor<dim>()
74+
:
75+
material_model_inputs.strain_rate[q]);
76+
77+
heating_model_outputs.heating_source_terms[q] = stress * deviatoric_strain_rate;
78+
79+
// If shear heating work fractions are provided, reduce the
80+
// overall heating by this amount (which is assumed to be converted into other forms of energy)
81+
if (shear_heating_out != nullptr)
82+
heating_model_outputs.heating_source_terms[q] *= shear_heating_out->shear_heating_work_fractions[q];
83+
84+
heating_model_outputs.lhs_latent_heat_terms[q] = 0.0;
85+
}
86+
}
87+
88+
89+
90+
template <int dim>
91+
void
92+
ShearHeatingAnisotropicViscosity<dim>::
93+
create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &material_model_outputs) const
94+
{
95+
const unsigned int n_points = material_model_outputs.viscosities.size();
96+
97+
if (material_model_outputs.template has_additional_output_object<MaterialModel::AnisotropicViscosity<dim>>() == false)
98+
{
99+
material_model_outputs.additional_outputs.push_back(
100+
std::make_unique<MaterialModel::AnisotropicViscosity<dim>> (n_points));
101+
}
102+
103+
this->get_material_model().create_additional_named_outputs(material_model_outputs);
104+
}
105+
}
106+
}
107+
108+
109+
110+
// explicit instantiations
111+
namespace aspect
112+
{
113+
namespace HeatingModel
114+
{
115+
ASPECT_REGISTER_HEATING_MODEL(ShearHeatingAnisotropicViscosity,
116+
"anisotropic shear heating",
117+
"Implementation of a standard model for shear heating. "
118+
"Adds the term: "
119+
"$ 2 \\eta \\left( \\varepsilon - \\frac{1}{3} \\text{tr} "
120+
"\\varepsilon \\mathbf 1 \\right) : \\left( \\varepsilon - \\frac{1}{3} "
121+
"\\text{tr} \\varepsilon \\mathbf 1 \\right)$ to the "
122+
"right-hand side of the temperature equation.")
123+
}
124+
}

0 commit comments

Comments
 (0)