Commit d9825ecb authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov

Merge branch 'TM_MPL' into 'master'

Use MPL in ThermoMechanics

See merge request ogs/ogs!3429
parents 1e33303e 2b799b14
......@@ -942,7 +942,7 @@ void ProjectData::parseProcesses(
name, *_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters,
_local_coordinate_system, integration_order,
process_config);
process_config, _media);
break;
case 3:
process = ProcessLib::ThermoMechanics::
......@@ -950,7 +950,7 @@ void ProjectData::parseProcesses(
name, *_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters,
_local_coordinate_system, integration_order,
process_config);
process_config, _media);
break;
}
}
......
......@@ -12,12 +12,14 @@
#include <cassert>
#include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
#include "MaterialLib/SolidModels/MechanicsBase.h"
#include "ParameterLib/Utils.h"
#include "ProcessLib/Output/CreateSecondaryVariables.h"
#include "ProcessLib/Utils/ProcessUtils.h"
#include "ThermoMechanicsProcess.h"
#include "ThermoMechanicsProcessData.h"
......@@ -25,17 +27,31 @@ namespace ProcessLib
{
namespace ThermoMechanics
{
void checkMPLProperties(
std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
{
std::array const required_solid_properties = {
MaterialPropertyLib::density, MaterialPropertyLib::thermal_expansivity,
MaterialPropertyLib::thermal_conductivity,
MaterialPropertyLib::specific_heat_capacity};
for (auto const& m : media)
{
checkRequiredProperties(m.second->phase("Solid"),
required_solid_properties);
}
}
template <int DisplacementDim>
std::unique_ptr<Process> createThermoMechanicsProcess(
std::string name,
MeshLib::Mesh& mesh,
std::string name, MeshLib::Mesh& mesh,
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<ProcessVariable> const& variables,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
boost::optional<ParameterLib::CoordinateSystem> const&
local_coordinate_system,
unsigned const integration_order,
BaseLib::ConfigTree const& config)
unsigned const integration_order, BaseLib::ConfigTree const& config,
std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
{
//! \ogs_file_param{prj__processes__process__type}
config.checkConfigParameter("type", "THERMO_MECHANICS");
......@@ -118,37 +134,6 @@ std::unique_ptr<Process> createThermoMechanicsProcess(
MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
parameters, local_coordinate_system, config);
// Reference solid density
auto const& reference_solid_density = ParameterLib::findParameter<double>(
config,
//! \ogs_file_param_special{prj__processes__process__THERMO_MECHANICS__reference_solid_density}
"reference_solid_density", parameters, 1, &mesh);
DBUG("Use '{:s}' as solid density parameter.",
reference_solid_density.name);
// Linear thermal expansion coefficient
auto const& linear_thermal_expansion_coefficient =
ParameterLib::findParameter<double>(
config,
//! \ogs_file_param_special{prj__processes__process__THERMO_MECHANICS__linear_thermal_expansion_coefficient}
"linear_thermal_expansion_coefficient", parameters, 1, &mesh);
DBUG("Use '{:s}' as linear thermal expansion coefficient.",
linear_thermal_expansion_coefficient.name);
// Specific heat capacity
auto const& specific_heat_capacity = ParameterLib::findParameter<double>(
config,
//! \ogs_file_param_special{prj__processes__process__THERMO_MECHANICS__specific_heat_capacity}
"specific_heat_capacity", parameters, 1, &mesh);
DBUG("Use '{:s}' as specific heat capacity parameter.",
specific_heat_capacity.name);
// Thermal conductivity // TODO To be changed as tensor input.
auto const& thermal_conductivity = ParameterLib::findParameter<double>(
config,
//! \ogs_file_param_special{prj__processes__process__THERMO_MECHANICS__thermal_conductivity}
"thermal_conductivity", parameters, 1, &mesh);
DBUG("Use '{:s}' as thermal conductivity parameter.",
thermal_conductivity.name);
// Specific body force
Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
{
......@@ -176,14 +161,17 @@ std::unique_ptr<Process> createThermoMechanicsProcess(
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value,
&mesh);
auto media_map =
MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
DBUG("Check the solid properties of ThermoMechanics process ...");
checkMPLProperties(media);
DBUG("Solid properties verified.");
ThermoMechanicsProcessData<DisplacementDim> process_data{
materialIDs(mesh),
std::move(media_map),
std::move(solid_constitutive_relations),
initial_stress,
reference_solid_density,
linear_thermal_expansion_coefficient,
specific_heat_capacity,
thermal_conductivity,
specific_body_force,
mechanics_process_id,
heat_conduction_process_id};
......@@ -208,7 +196,8 @@ template std::unique_ptr<Process> createThermoMechanicsProcess<2>(
boost::optional<ParameterLib::CoordinateSystem> const&
local_coordinate_system,
unsigned const integration_order,
BaseLib::ConfigTree const& config);
BaseLib::ConfigTree const& config,
std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
template std::unique_ptr<Process> createThermoMechanicsProcess<3>(
std::string name,
......@@ -219,7 +208,8 @@ template std::unique_ptr<Process> createThermoMechanicsProcess<3>(
boost::optional<ParameterLib::CoordinateSystem> const&
local_coordinate_system,
unsigned const integration_order,
BaseLib::ConfigTree const& config);
BaseLib::ConfigTree const& config,
std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
} // namespace ThermoMechanics
} // namespace ProcessLib
......@@ -11,6 +11,7 @@
#pragma once
#include <boost/optional.hpp>
#include <map>
#include <memory>
#include <vector>
......@@ -22,6 +23,12 @@ namespace MeshLib
{
class Mesh;
}
namespace MaterialPropertyLib
{
class Medium;
}
namespace ParameterLib
{
struct CoordinateSystem;
......@@ -40,37 +47,34 @@ namespace ThermoMechanics
{
template <int DisplacementDim>
std::unique_ptr<Process> createThermoMechanicsProcess(
std::string name,
MeshLib::Mesh& mesh,
std::string name, MeshLib::Mesh& mesh,
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<ProcessVariable> const& variables,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
boost::optional<ParameterLib::CoordinateSystem> const&
local_coordinate_system,
unsigned const integration_order,
BaseLib::ConfigTree const& config);
unsigned const integration_order, BaseLib::ConfigTree const& config,
std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
extern template std::unique_ptr<Process> createThermoMechanicsProcess<2>(
std::string name,
MeshLib::Mesh& mesh,
std::string name, MeshLib::Mesh& mesh,
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<ProcessVariable> const& variables,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
boost::optional<ParameterLib::CoordinateSystem> const&
local_coordinate_system,
unsigned const integration_order,
BaseLib::ConfigTree const& config);
unsigned const integration_order, BaseLib::ConfigTree const& config,
std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
extern template std::unique_ptr<Process> createThermoMechanicsProcess<3>(
std::string name,
MeshLib::Mesh& mesh,
std::string name, MeshLib::Mesh& mesh,
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<ProcessVariable> const& variables,
std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
boost::optional<ParameterLib::CoordinateSystem> const&
local_coordinate_system,
unsigned const integration_order,
BaseLib::ConfigTree const& config);
unsigned const integration_order, BaseLib::ConfigTree const& config,
std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
} // namespace ThermoMechanics
} // namespace ProcessLib
......@@ -71,7 +71,7 @@ AddTest(
RUNTIME 17
DIFF_DATA
tm1_2Dbeam_ts_1_t_1.000000.vtu tm1_2Dbeam_ts_1_t_1.000000.vtu temperature temperature 1e-11 0.0
tm1_2Dbeam_ts_1_t_1.000000.vtu tm1_2Dbeam_ts_1_t_1.000000.vtu displacement displacement 1e-11 0.0
tm1_2Dbeam_ts_1_t_1.000000.vtu tm1_2Dbeam_ts_1_t_1.000000.vtu displacement displacement 5e-11 0.0
)
AddTest(
......@@ -113,7 +113,7 @@ AddTest(
RUNTIME 17
DIFF_DATA
tm1_3Dgravity_ts_1_t_1.000000.vtu tm1_3Dgravity_ts_1_t_1.000000.vtu temperature temperature 1e-11 0.0
tm1_3Dgravity_ts_1_t_1.000000.vtu tm1_3Dgravity_ts_1_t_1.000000.vtu displacement displacement 1e-11 0.0
tm1_3Dgravity_ts_1_t_1.000000.vtu tm1_3Dgravity_ts_1_t_1.000000.vtu displacement displacement 1e-7 0.0
)
AddTest(
......@@ -155,9 +155,9 @@ AddTest(
REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
RUNTIME 17
DIFF_DATA
tm2_1D1bt_ts_50_t_5.000000.vtu tm2_1D1bt_ts_50_t_5.000000.vtu temperature temperature 1e-10 0.0
tm2_1D1bt_ts_50_t_5.000000.vtu tm2_1D1bt_ts_50_t_5.000000.vtu temperature temperature 5e-8 0.0
tm2_1D1bt_ts_50_t_5.000000.vtu tm2_1D1bt_ts_50_t_5.000000.vtu displacement displacement 1e-10 0.0
tm2_1D1bt_ts_100_t_10.000000.vtu tm2_1D1bt_ts_100_t_10.000000.vtu temperature temperature 1e-10 0.0
tm2_1D1bt_ts_100_t_10.000000.vtu tm2_1D1bt_ts_100_t_10.000000.vtu temperature temperature 3e-7 0.0
tm2_1D1bt_ts_100_t_10.000000.vtu tm2_1D1bt_ts_100_t_10.000000.vtu displacement displacement 1e-10 0.0
)
......@@ -251,10 +251,10 @@ AddTest(
TESTER vtkdiff
REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
DIFF_DATA
expected_tm_a_quad_ts_20_t_20000.000000.vtu tm_a_quad_ts_20_t_20000.000000.vtu displacement displacement 5e-10 1e-15
expected_tm_a_quad_ts_20_t_20000.000000.vtu tm_a_quad_ts_20_t_20000.000000.vtu temperature temperature 5e-10 1e-7
expected_tm_a_quad_ts_20_t_20000.000000.vtu tm_a_quad_ts_20_t_20000.000000.vtu sigma sigma 1e-6 0
expected_tm_a_quad_ts_20_t_20000.000000.vtu tm_a_quad_ts_20_t_20000.000000.vtu epsilon epsilon 5e-10 1e-8
expected_tm_a_quad_ts_20_t_20000.000000.vtu tm_a_quad_ts_20_t_20000.000000.vtu displacement displacement 2e-7 1e-15
expected_tm_a_quad_ts_20_t_20000.000000.vtu tm_a_quad_ts_20_t_20000.000000.vtu temperature temperature 5e-4 1e-7
expected_tm_a_quad_ts_20_t_20000.000000.vtu tm_a_quad_ts_20_t_20000.000000.vtu sigma sigma 5e-4 0
expected_tm_a_quad_ts_20_t_20000.000000.vtu tm_a_quad_ts_20_t_20000.000000.vtu epsilon epsilon 1e-8 1e-8
)
AddTest(
......@@ -294,9 +294,9 @@ AddTest(
REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
RUNTIME 24
DIFF_DATA
ExpectedCreepAfterExcavation_ts_61_t_4320000.000000.vtu CreepAfterExcavation_ts_61_t_4320000.000000.vtu sigma sigma 5e-6 0
ExpectedCreepAfterExcavation_ts_61_t_4320000.000000.vtu CreepAfterExcavation_ts_61_t_4320000.000000.vtu epsilon epsilon 1e-15 0
ExpectedCreepAfterExcavation_ts_61_t_4320000.000000.vtu CreepAfterExcavation_ts_61_t_4320000.000000.vtu displacement displacement 1e-16 1e-9
CreepAfterExcavation_ts_61_t_4320000.000000.vtu CreepAfterExcavation_ts_61_t_4320000.000000.vtu sigma sigma 5e-6 0
CreepAfterExcavation_ts_61_t_4320000.000000.vtu CreepAfterExcavation_ts_61_t_4320000.000000.vtu epsilon epsilon 1e-15 0
CreepAfterExcavation_ts_61_t_4320000.000000.vtu CreepAfterExcavation_ts_61_t_4320000.000000.vtu displacement displacement 1e-16 1e-9
)
# Basic test that MFront models work for TM.
......
......@@ -12,8 +12,11 @@
#pragma once
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/MPL/Property.h"
#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
#include "MaterialLib/MPL/Utils/FormKelvinVectorFromThermalExpansivity.h"
#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
#include "ThermoMechanicsFEM.h"
namespace ProcessLib
{
......@@ -71,9 +74,6 @@ ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
ip_data.eps_m_prev.setZero(kelvin_vector_size);
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
ip_data.solid_density =
_process_data.reference_solid_density(0, x_position)[0];
ip_data.solid_density_prev = ip_data.solid_density;
ip_data.N = shape_matrices[ip].N;
ip_data.dNdx = shape_matrices[ip].dNdx;
......@@ -168,6 +168,9 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
auto const& medium = _process_data.media_map->getMedium(_element.getID());
auto const& solid_phase = medium->phase("Solid");
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
......@@ -195,27 +198,27 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
auto& state = _ip_data[ip].material_state_variables;
const double T_ip = N.dot(T); // T at integration point
double const dT = N.dot(T_dot) * dt;
// calculate thermally induced strain
// assume isotropic thermal expansion
auto const alpha = _process_data.linear_thermal_expansion_coefficient(
t, x_position)[0];
double const linear_thermal_strain_increment = alpha * dT;
// Consider also anisotropic thermal expansion.
auto const solid_linear_thermal_expansivity_vector =
MPL::formKelvinVectorFromThermalExpansivity<DisplacementDim>(
solid_phase
.property(
MaterialPropertyLib::PropertyType::thermal_expansivity)
.value(variables, x_position, t, dt));
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
dthermal_strain =
solid_linear_thermal_expansivity_vector.eval() * dT;
//
// displacement equation, displacement part
//
eps.noalias() = B * u;
using Invariants = MathLib::KelvinVector::Invariants<
MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value>;
// assume isotropic thermal expansion
const double T_ip = N.dot(T); // T at integration point
eps_m.noalias() =
eps_m_prev + eps - eps_prev -
linear_thermal_strain_increment * Invariants::identity2;
eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
variables_prev[static_cast<int>(MPL::Variable::stress)]
.emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
......@@ -260,14 +263,9 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
.noalias() = N;
}
// calculate real density
// rho_s_{n+1} * (V_{n} + dV) = rho_s_n * V_n
// dV = 3 * alpha * dT * V_0
// rho_s_{n+1} = rho_s_n / (1 + 3 * alpha * dT )
// see reference solid density description for details.
auto& rho_s = _ip_data[ip].solid_density;
rho_s = _ip_data[ip].solid_density_prev /
(1 + 3 * linear_thermal_strain_increment);
auto const rho_s =
solid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
auto const& b = _process_data.specific_body_force;
local_rhs.template block<displacement_size, 1>(displacement_index, 0)
......@@ -276,9 +274,12 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
//
// displacement equation, temperature part
//
KuT.noalias() +=
B.transpose() * C * alpha * Invariants::identity2 * N * w;
// The computation of KuT can be ignored.
auto const alpha_T_tensor =
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
solid_linear_thermal_expansivity_vector.eval());
KuT.noalias() += B.transpose() * (C * alpha_T_tensor.eval()) * N * w;
if (_ip_data[ip].solid_material.getConstitutiveModel() ==
MaterialLib::Solids::ConstitutiveModel::CreepBGRa)
{
......@@ -294,10 +295,21 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
// temperature equation, temperature part;
//
auto const lambda =
_process_data.thermal_conductivity(t, x_position)[0];
KTT.noalias() += dNdx.transpose() * lambda * dNdx * w;
solid_phase
.property(
MaterialPropertyLib::PropertyType::thermal_conductivity)
.value(variables, x_position, t, dt);
GlobalDimMatrixType const thermal_conductivity =
MaterialPropertyLib::formEigenTensor<DisplacementDim>(lambda);
auto const c = _process_data.specific_heat_capacity(t, x_position)[0];
KTT.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
auto const c =
solid_phase
.property(
MaterialPropertyLib::PropertyType::specific_heat_capacity)
.template value<double>(variables, x_position, t, dt);
DTT.noalias() += N.transpose() * rho_s * c * N * w;
}
......@@ -376,6 +388,8 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
MPL::VariableArray variables_prev;
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
auto const& medium = _process_data.media_map->getMedium(_element.getID());
auto const& solid_phase = medium->phase("Solid");
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
......@@ -406,25 +420,27 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
const double T_ip = N.dot(local_T); // T at integration point
double const dT_ip = N.dot(local_Tdot) * dt;
// calculate thermally induced strain
// assume isotropic thermal expansion
auto const alpha = _process_data.linear_thermal_expansion_coefficient(
t, x_position)[0];
double const linear_thermal_strain_increment = alpha * dT_ip;
variables[static_cast<int>(MPL::Variable::temperature)].emplace<double>(
T_ip);
//
// displacement equation, displacement part
//
eps.noalias() = B * local_u;
using Invariants = MathLib::KelvinVector::Invariants<
MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value>;
// Consider also anisotropic thermal expansion.
auto const solid_linear_thermal_expansivity_vector =
MPL::formKelvinVectorFromThermalExpansivity<DisplacementDim>(
solid_phase
.property(
MaterialPropertyLib::PropertyType::thermal_expansivity)
.value(variables, x_position, t, dt));
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
dthermal_strain =
solid_linear_thermal_expansivity_vector.eval() * dT_ip;
// assume isotropic thermal expansion
eps_m.noalias() =
eps_m_prev + eps - eps_prev -
linear_thermal_strain_increment * Invariants::identity2;
eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
variables_prev[static_cast<int>(MPL::Variable::stress)]
.emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
......@@ -437,8 +453,6 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
variables[static_cast<int>(MPL::Variable::mechanical_strain)]
.emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
eps_m);
variables[static_cast<int>(MPL::Variable::temperature)].emplace<double>(
T_ip);
auto&& solution = _ip_data[ip].solid_material.integrateStress(
variables_prev, variables, t, x_position, dt, *state);
......@@ -466,14 +480,9 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
.noalias() = N;
}
// calculate real density
// rho_s_{n+1} * (V_{n} + dV) = rho_s_n * V_n
// dV = 3 * alpha * dT * V_0
// rho_s_{n+1} = rho_s_n / (1 + 3 * alpha * dT )
// see reference solid density description for details.
auto& rho_s = _ip_data[ip].solid_density;
rho_s = _ip_data[ip].solid_density_prev /
(1 + 3 * linear_thermal_strain_increment);
auto const rho_s =
solid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
auto const& b = _process_data.specific_body_force;
local_rhs.noalias() -=
......@@ -516,6 +525,9 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
auto const& medium = _process_data.media_map->getMedium(_element.getID());
auto const& solid_phase = medium->phase("Solid");
MPL::VariableArray variables;
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
......@@ -524,27 +536,29 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
auto const& N = _ip_data[ip].N;
auto const& dNdx = _ip_data[ip].dNdx;
// calculate real density
// rho_s_{n+1} * (V_{n} + dV) = rho_s_n * V_n
// dV = 3 * alpha * dT * V_0
// rho_s_{n+1} = rho_s_n / (1 + 3 * alpha * dT )
// see reference solid density description for details.
auto& rho_s = _ip_data[ip].solid_density;
// calculate thermally induced strain
// assume isotropic thermal expansion
auto const alpha = _process_data.linear_thermal_expansion_coefficient(
t, x_position)[0];
double const dT_ip = N.dot(local_dT);
double const linear_thermal_strain_increment = alpha * dT_ip;
rho_s = _ip_data[ip].solid_density_prev /
(1 + 3 * linear_thermal_strain_increment);
auto const c_p = _process_data.specific_heat_capacity(t, x_position)[0];
const double T_ip = N.dot(local_T); // T at integration point
variables[static_cast<int>(MPL::Variable::temperature)].emplace<double>(
T_ip);
auto const rho_s =
solid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
auto const c_p =
solid_phase.property(MPL::PropertyType::specific_heat_capacity)
.template value<double>(variables, x_position, t, dt);
mass.noalias() += N.transpose() * rho_s * c_p * N * w;
auto const lambda =
_process_data.thermal_conductivity(t, x_position)[0];
laplace.noalias() += dNdx.transpose() * lambda * dNdx * w;
solid_phase
.property(
MaterialPropertyLib::PropertyType::thermal_conductivity)
.value(variables, x_position, t, dt);
GlobalDimMatrixType const thermal_conductivity =
MaterialPropertyLib::formEigenTensor<DisplacementDim>(lambda);
laplace.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
}
local_Jac.noalias() += laplace + mass / dt;
......
......@@ -58,8 +58,6 @@ struct IntegrationPointData final
std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
DisplacementDim>::MaterialStateVariables>
material_state_variables;
double solid_density;
double solid_density_prev;
double integration_weight;
typename ShapeMatricesType::NodalRowVectorType N;
......@@ -70,7 +68,6 @@ struct IntegrationPointData final
eps_prev = eps;
eps_m_prev = eps_m;
sigma_prev = sigma;
solid_density_prev = solid_density;
material_state_variables->pushBackState();
}
......@@ -85,6 +82,46 @@ struct SecondaryData
std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
};
/**
* \brief Local assembler of ThermoMechanics process.
*
* In this assembler, the solid density can be expressed as an exponential
* function of temperature, if a temperature dependent solid density is
* assumed.
* The theory behind this is given below.
*
* During the temperature variation, the solid mass balance is
* \f[
* \rho_s^{n+1} V^{n+1} = \rho_s^n V^n,
* \f]
* with \f$\rho_s\f$ the density, \f$V\f$, the volume, and \f$n\f$ the index
* of the time step.
* Under pure thermo-mechanics condition, the volume change along with
* the temperature change is given by
* \f[
* V^{n+1} = V^{n} + \alpha_T dT V^{n} = (1+\alpha_T dT)V^{n},
* \f]
* where \f$\alpha_T\f$ is the volumetric thermal expansivity of the solid.