From e48ccd60d9adde40c5b5b9b2cb6ede1b5a91987f Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Mon, 26 Aug 2019 17:06:27 +0200
Subject: [PATCH] [PL] TM; Add initial stress evaluation.
---
.../THERMO_MECHANICS/t_initial_stress.md | 1 +
.../CreateThermoMechanicsProcess.cpp | 33 ++++++++++++++--
.../ThermoMechanics/ThermoMechanicsFEM-impl.h | 39 +++++--------------
.../ThermoMechanics/ThermoMechanicsFEM.h | 22 ++++++++++-
.../ThermoMechanicsProcessData.h | 4 ++
.../non_equilibrium_initial_state.prj | 5 +--
6 files changed, 65 insertions(+), 39 deletions(-)
create mode 100644 Documentation/ProjectFile/prj/processes/process/THERMO_MECHANICS/t_initial_stress.md
diff --git a/Documentation/ProjectFile/prj/processes/process/THERMO_MECHANICS/t_initial_stress.md b/Documentation/ProjectFile/prj/processes/process/THERMO_MECHANICS/t_initial_stress.md
new file mode 100644
index 00000000000..91a430a2e69
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/THERMO_MECHANICS/t_initial_stress.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::ThermoMechanics::ThermoMechanicsProcessData::initial_stress
diff --git a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
index 7b30e5574a9..deaf6f9c047 100644
--- a/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
+++ b/ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
@@ -167,11 +167,36 @@ std::unique_ptr<Process> createThermoMechanicsProcess(
std::copy_n(b.data(), b.size(), specific_body_force.data());
}
+ // Initial stress conditions
+ ParameterLib::Parameter<double> const* initial_stress = nullptr;
+ auto const initial_stress_parameter_name =
+ //! \ogs_file_param{prj__processes__process__THERMO_MECHANICS__initial_stress}
+ config.getConfigParameterOptional<std::string>("initial_stress");
+ if (initial_stress_parameter_name)
+ {
+ initial_stress = ParameterLib::findParameterOptional<double>(
+ *initial_stress_parameter_name, parameters,
+ // Symmetric tensor size, 4 or 6, not a Kelvin vector.
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value,
+ &mesh);
+ }
+ if (initial_stress != nullptr)
+ {
+ DBUG("Use '%s' as initial stress parameter.",
+ initial_stress->name.c_str());
+ }
+
ThermoMechanicsProcessData<DisplacementDim> process_data{
- materialIDs(mesh), std::move(solid_constitutive_relations),
- reference_solid_density, linear_thermal_expansion_coefficient,
- specific_heat_capacity, thermal_conductivity,
- specific_body_force, mechanics_process_id,
+ materialIDs(mesh),
+ 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};
SecondaryVariableCollection secondary_variables;
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
index 18456690f5d..9e472c3a692 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
@@ -59,10 +59,14 @@ ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
static const int kelvin_vector_size =
MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value;
+ // Initialize current time step values
ip_data.sigma.setZero(kelvin_vector_size);
- ip_data.sigma_prev.setZero(kelvin_vector_size);
ip_data.eps.setZero(kelvin_vector_size);
- ip_data.eps_prev.setZero(kelvin_vector_size);
+
+ // Previous time step values are not initialized and are set later.
+ ip_data.sigma_prev.resize(kelvin_vector_size);
+ ip_data.eps_prev.resize(kelvin_vector_size);
+
ip_data.eps_m.setZero(kelvin_vector_size);
ip_data.eps_m_prev.setZero(kelvin_vector_size);
ParameterLib::SpatialPosition x_position;
@@ -74,28 +78,6 @@ ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
ip_data.dNdx = shape_matrices[ip].dNdx;
_secondary_data.N[ip] = shape_matrices[ip].N;
-
- // Initialize current time step values
- if (_process_data.nonequilibrium_stress)
- {
- // Computation of non-equilibrium stress.
- x_position.setCoordinates(MathLib::Point3d(
- interpolateCoordinates<ShapeFunction, ShapeMatricesType>(
- e, ip_data.N)));
- std::vector<double> sigma_neq_data = (*_process_data
- .nonequilibrium_stress)(
- std::numeric_limits<double>::quiet_NaN() /* time independent */,
- x_position);
- ip_data.sigma_neq =
- Eigen::Map<typename BMatricesType::KelvinVectorType>(
- sigma_neq_data.data(),
- MathLib::KelvinVector::KelvinVectorDimensions<
- DisplacementDim>::value,
- 1);
- }
- // Initialization from non-equilibrium sigma, which is zero by
- // default, or is set to some value.
- ip_data.sigma = ip_data.sigma_neq;
}
}
@@ -203,7 +185,6 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
auto& sigma = _ip_data[ip].sigma;
auto const& sigma_prev = _ip_data[ip].sigma_prev;
- auto const& sigma_neq = _ip_data[ip].sigma_neq;
auto& eps = _ip_data[ip].eps;
auto const& eps_prev = _ip_data[ip].eps_prev;
@@ -275,9 +256,8 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
auto const& b = _process_data.specific_body_force;
local_rhs.template block<displacement_size, 1>(displacement_index, 0)
- .noalias() -= (B.transpose() * (sigma - sigma_neq) -
- N_u.transpose() * rho_s * b) *
- w;
+ .noalias() -=
+ (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
//
// displacement equation, temperature part
@@ -287,8 +267,7 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
if (_ip_data[ip].solid_material.getConstitutiveModel() ==
MaterialLib::Solids::ConstitutiveModel::CreepBGRa)
{
- auto const s =
- Invariants::deviatoric_projection * (sigma - sigma_neq);
+ auto const s = Invariants::deviatoric_projection * sigma;
double const norm_s = Invariants::FrobeniusNorm(s);
const double creep_coefficient =
_ip_data[ip].solid_material.getTemperatureRelatedCoefficient(
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
index 3bfa7c764fe..fa7a10e66ca 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
@@ -151,10 +151,28 @@ public:
{
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
-
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- _ip_data[ip].pushBackState();
+ auto& ip_data = _ip_data[ip];
+
+ /// Set initial stress from parameter.
+ if (_process_data.initial_stress != nullptr)
+ {
+ ParameterLib::SpatialPosition const x_position{
+ boost::none, _element.getID(), ip,
+ MathLib::Point3d(interpolateCoordinates<ShapeFunction,
+ ShapeMatricesType>(
+ _element, ip_data.N))};
+
+ ip_data.sigma =
+ MathLib::KelvinVector::symmetricTensorToKelvinVector<
+ DisplacementDim>((*_process_data.initial_stress)(
+ std::numeric_limits<
+ double>::quiet_NaN() /* time independent */,
+ x_position));
+ }
+
+ ip_data.pushBackState();
}
}
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
index eee568d37a9..7ee8e144415 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
@@ -38,6 +38,10 @@ struct ThermoMechanicsProcessData
MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
solid_materials;
+ /// Optional, initial stress field. A symmetric tensor, short vector
+ /// representation of length 4 or 6, ParameterLib::Parameter<double>.
+ ParameterLib::Parameter<double> const* const initial_stress;
+
ParameterLib::Parameter<double> const& reference_solid_density;
ParameterLib::Parameter<double> const& linear_thermal_expansion_coefficient;
ParameterLib::Parameter<double> const& specific_heat_capacity;
diff --git a/Tests/Data/ThermoMechanics/InitialStates/non_equilibrium_initial_state.prj b/Tests/Data/ThermoMechanics/InitialStates/non_equilibrium_initial_state.prj
index 4e3b8bfe405..79bdbc0c4ae 100644
--- a/Tests/Data/ThermoMechanics/InitialStates/non_equilibrium_initial_state.prj
+++ b/Tests/Data/ThermoMechanics/InitialStates/non_equilibrium_initial_state.prj
@@ -30,15 +30,14 @@
<secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
<secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
</secondary_variables>
- <nonequilibrium_state_variables>
- <stress>nonequilibrium_stress</stress>
- </nonequilibrium_state_variables>
+ <initial_stress>nonequilibrium_stress</initial_stress>
</process>
</processes>
<time_loop>
<processes>
<process ref="TM">
<nonlinear_solver>basic_newton</nonlinear_solver>
+ <compensate_initial_out_of_balance_forces>true</compensate_initial_out_of_balance_forces>
<convergence_criterion>
<type>DeltaX</type>
<norm_type>NORM2</norm_type>
--
GitLab