diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index d55b03082197bd04cb94d4f28660c82114a8923b..0965f8009c6c145b45c989ca9f5f97f00d7c1fed 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -137,6 +137,7 @@ struct ConstitutiveTempData
ViscosityData viscosity_data;
PhaseTransitionData phase_transition_data;
PorosityDerivativeData porosity_d_data;
+ SolidDensityDerivativeData solid_density_d_data;
using DisplacementDimVector = Eigen::Matrix<double, DisplacementDim, 1>;
using DisplacementDimMatrix =
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index bb9847ebdc66c654b777f5abc725fdddb1b6db26..b9526cdf6d34197e8049e328e2930dc1ec84772c 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -19,6 +19,7 @@
#include "PureLiquidDensity.h"
#include "Saturation.h"
#include "SolidCompressibility.h"
+#include "SolidDensity.h"
#include "SolidMechanics.h"
#include "SolidThermalExpansion.h"
#include "Swelling.h"
@@ -61,10 +62,13 @@ struct ConstitutiveModels
PhaseTransitionModel const& phase_transition_model;
#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
PorosityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>
+ porosity_model;
+ SolidDensityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>
+ solid_density_model;
#else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- PorosityModel
+ PorosityModel porosity_model;
+ SolidDensityModel solid_density_model;
#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- porosity_model;
};
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9978c1516908253b798ee355920cc70fb4c5809b
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.cpp
@@ -0,0 +1,79 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "SolidDensity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+
+void SolidDensityModel::eval(
+ SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ TemperatureData const& T_data,
+ SolidDensityData& solid_density_data,
+ SolidDensityDerivativeData& solid_density_d_data) const
+{
+ MaterialPropertyLib::VariableArray variables;
+ variables.temperature = T_data.T;
+
+ auto const& mpl_solid_density =
+ media_data.solid[MaterialPropertyLib::PropertyType::density];
+
+ solid_density_data.rho_SR = mpl_solid_density.template value<double>(
+ variables, x_t.x, x_t.t, x_t.dt);
+
+ solid_density_d_data.drho_SR_dT = mpl_solid_density.template dValue<double>(
+ variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t,
+ x_t.dt);
+}
+
+template <int DisplacementDim>
+void SolidDensityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>::
+ eval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ TemperatureData const& T_data,
+ BiotData const& biot,
+ StrainData<DisplacementDim> const& strain_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ SolidDensityData& solid_density_data,
+ SolidDensityDerivativeData& solid_density_d_data) const
+{
+ MaterialPropertyLib::VariableArray variables;
+ variables.temperature = T_data.T;
+
+ static int const KelvinVectorSize =
+ MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+ using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+ double const div_u = Invariants::trace(strain_data.eps);
+
+ auto const& mpl_solid_density =
+ media_data.solid[MaterialPropertyLib::PropertyType::density];
+
+ auto const rho_ref_SR = mpl_solid_density.template value<double>(
+ variables, x_t.x, x_t.t, x_t.dt);
+
+ solid_density_data.rho_SR =
+ rho_ref_SR *
+ (1. - s_therm_exp_data.thermal_volume_strain + (biot() - 1.) * div_u);
+
+ solid_density_d_data.drho_SR_dT =
+ mpl_solid_density.template dValue<double>(
+ variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t,
+ x_t.dt) *
+ (1. - s_therm_exp_data.thermal_volume_strain +
+ (biot() - 1.) * div_u) -
+ rho_ref_SR * s_therm_exp_data.beta_T_SR;
+}
+
+template struct SolidDensityModelNonConstantSolidPhaseVolumeFraction<2>;
+template struct SolidDensityModelNonConstantSolidPhaseVolumeFraction<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h
index 399f73851f0158ac2155ede0bd0ef5fc03ddc5a8..ad14ab20ff85c3de00210a6740324b59e3866e4f 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h
@@ -10,12 +10,20 @@
#pragma once
#include "Base.h"
+#include "Biot.h"
+#include "ProcessLib/ConstitutiveRelations/StrainData.h"
#include "ProcessLib/Reflection/ReflectionData.h"
+#include "SolidThermalExpansion.h"
namespace ProcessLib::TH2M
{
namespace ConstitutiveRelations
{
+struct SolidDensityDerivativeData
+{
+ double drho_SR_dT = nan;
+};
+
struct SolidDensityData
{
double rho_SR = nan;
@@ -29,5 +37,31 @@ struct SolidDensityData
R::makeReflectionData("solid_density", &Self::rho_SR)};
}
};
+
+struct SolidDensityModel
+{
+ void eval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ TemperatureData const& T_data,
+ SolidDensityData& solid_density_data,
+ SolidDensityDerivativeData& solid_density_d_data) const;
+};
+
+template <int DisplacementDim>
+struct SolidDensityModelNonConstantSolidPhaseVolumeFraction
+{
+ void eval(
+ SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ TemperatureData const& T_data,
+ BiotData const& biot,
+ StrainData<DisplacementDim> const& strain_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ SolidDensityData& solid_density_data,
+ SolidDensityDerivativeData& solid_density_d_data) const;
+};
+
+extern template struct SolidDensityModelNonConstantSolidPhaseVolumeFraction<2>;
+extern template struct SolidDensityModelNonConstantSolidPhaseVolumeFraction<3>;
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index b4a80dc60bd3bd8985c07e63f5d8489f18230553..341e92cbee58aacf8e8973a26fe3dd294a233977 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -239,6 +239,13 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
ip_out.porosity_data, ip_cv.porosity_d_data);
+ models.solid_density_model.eval(
+ {pos, t, dt}, media_data, T_data,
+#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+ ip_cv.biot_data, ip_out.eps_data, ip_cv.s_therm_exp_data,
+#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+ ip_out.solid_density_data, ip_cv.solid_density_d_data);
+
MPL::VariableArray vars;
MPL::VariableArray vars_prev;
vars.temperature = T;
@@ -254,20 +261,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
vars.porosity = ip_out.porosity_data.phi;
- // solid phase density
- auto const rho_ref_SR =
- solid_phase.property(MPL::PropertyType::density)
- .template value<double>(
- vars, pos, t, std::numeric_limits<double>::quiet_NaN());
-
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- ip_out.solid_density_data.rho_SR =
- rho_ref_SR * (1. - ip_cv.s_therm_exp_data.thermal_volume_strain +
- (ip_cv.biot_data() - 1.) * div_u);
-#else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- ip_out.solid_density_data.rho_SR = rho_ref_SR;
-#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
-
auto const& c = ip_cv.phase_transition_data;
auto const phi_L =
@@ -348,23 +341,13 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
liquid_phase.property(MPL::PropertyType::density)
.template dValue<double>(vars, MPL::Variable::temperature, pos,
t, dt);
- auto const drho_SR_dT =
- solid_phase.property(MPL::PropertyType::density)
- .template dValue<double>(vars, MPL::Variable::temperature,
- pos, t, dt)
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- * (1. - ip_cv.s_therm_exp_data.thermal_volume_strain +
- (ip_cv.biot_data() - 1.) * div_u) -
- rho_ref_SR * ip_cv.s_therm_exp_data.beta_T_SR
-#endif
- ;
ip_cv.drho_u_eff_dT =
phi_G * c.drho_GR_dT * c.uG +
phi_G * ip_out.fluid_density_data.rho_GR * c.du_G_dT +
phi_L * drho_LR_dT * c.uL +
phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dT +
- phi_S * drho_SR_dT * u_S +
+ phi_S * ip_cv.solid_density_d_data.drho_SR_dT * u_S +
phi_S * ip_out.solid_density_data.rho_SR * cpS +
ip_cv.porosity_d_data.dphi_S_dT * ip_out.solid_density_data.rho_SR *
u_S;
@@ -469,7 +452,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
phi_L * ip_out.fluid_density_data.rho_LR * c.dh_L_dT +
ip_cv.porosity_d_data.dphi_S_dT * ip_out.solid_density_data.rho_SR *
ip_data.h_S +
- phi_S * drho_SR_dT * ip_data.h_S +
+ phi_S * ip_cv.solid_density_d_data.drho_SR_dT * ip_data.h_S +
phi_S * ip_out.solid_density_data.rho_SR * cpS;
ip_cv.drho_u_eff_dp_GR =