diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Base.h b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
index 7ee7b82d247822c6e4579e2bd4988fc797332889..d0f26028eca002d68e9daa55417801abd26c1742 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Base.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
@@ -51,6 +51,8 @@ struct TemperatureData
double T_prev;
};
+using ReferenceTemperatureData =
+ BaseLib::StrongType<double, struct ReferenceTemperatureTag>;
using GasPressureData = BaseLib::StrongType<double, struct GasPressureTag>;
using CapillaryPressureData =
BaseLib::StrongType<double, struct CapillaryPressureTag>;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index a25ddf0bf9c693fd34eb651da0dd14d1ce021a89..0965f8009c6c145b45c989ca9f5f97f00d7c1fed 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -20,12 +20,14 @@
#include "MechanicalStrain.h"
#include "PermeabilityData.h"
#include "PhaseTransitionData.h"
+#include "Porosity.h"
#include "ProcessLib/ConstitutiveRelations/StrainData.h"
#include "ProcessLib/ConstitutiveRelations/StressData.h"
#include "ProcessLib/Reflection/ReflectionData.h"
#include "PureLiquidDensity.h"
#include "Saturation.h"
#include "SolidCompressibility.h"
+#include "SolidDensity.h"
#include "SolidMechanics.h"
#include "SolidThermalExpansion.h"
#include "Swelling.h"
@@ -93,6 +95,8 @@ struct OutputData
MassMoleFractionsData mass_mole_fractions_data;
FluidDensityData fluid_density_data;
VapourPartialPressureData vapour_pressure_data;
+ PorosityData porosity_data;
+ SolidDensityData solid_density_data;
static auto reflect()
{
@@ -103,7 +107,9 @@ struct OutputData
&Self::enthalpy_data,
&Self::mass_mole_fractions_data,
&Self::fluid_density_data,
- &Self::vapour_pressure_data);
+ &Self::vapour_pressure_data,
+ &Self::porosity_data,
+ &Self::solid_density_data);
}
};
@@ -130,6 +136,8 @@ struct ConstitutiveTempData
EquivalentPlasticStrainData equivalent_plastic_strain_data;
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 7fe33f3683b5b005ee2e94dd4cb2eb120f495ae2..b9526cdf6d34197e8049e328e2930dc1ec84772c 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -15,9 +15,11 @@
#include "MechanicalStrain.h"
#include "PermeabilityModel.h"
#include "PhaseTransitionModel.h"
+#include "Porosity.h"
#include "PureLiquidDensity.h"
#include "Saturation.h"
#include "SolidCompressibility.h"
+#include "SolidDensity.h"
#include "SolidMechanics.h"
#include "SolidThermalExpansion.h"
#include "Swelling.h"
@@ -58,6 +60,15 @@ struct ConstitutiveModels
PermeabilityModel<DisplacementDim> permeability_model;
PureLiquidDensityModel pure_liquid_density_model;
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 porosity_model;
+ SolidDensityModel solid_density_model;
+#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
};
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..508dc689d3a8c32a338a0cb77093b0a259c75a66
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp
@@ -0,0 +1,77 @@
+/**
+ * \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 "Porosity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+
+void PorosityModel::eval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ PorosityData& porosity_data,
+ PorosityDerivativeData& porosity_d_data) const
+{
+ MaterialPropertyLib::VariableArray variables;
+
+ auto const& mpl_porosity =
+ media_data.medium[MaterialPropertyLib::PropertyType::porosity];
+
+ porosity_data.phi =
+ mpl_porosity.template value<double>(variables, x_t.x, x_t.t, x_t.dt);
+
+ porosity_d_data.dphi_dT = mpl_porosity.template dValue<double>(
+ variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t,
+ x_t.dt);
+}
+
+template <int DisplacementDim>
+void PorosityModelNonConstantSolidPhaseVolumeFraction<DisplacementDim>::eval(
+ SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ BiotData const& biot,
+ StrainData<DisplacementDim> const& strain_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ PorosityData& porosity_data,
+ PorosityDerivativeData& porosity_d_data) const
+{
+ MaterialPropertyLib::VariableArray variables;
+
+ auto const& mpl_porosity =
+ media_data.medium[MaterialPropertyLib::PropertyType::porosity];
+
+ double const phi_0 =
+ mpl_porosity.template value<double>(variables, x_t.x, x_t.t, x_t.dt);
+
+ 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);
+
+ double const phi_S =
+ (1. - phi_0) *
+ (1. + s_therm_exp_data.thermal_volume_strain - biot() * div_u);
+
+ porosity_data.phi = 1. - phi_S;
+
+ auto const dphi_0_dT = mpl_porosity.template dValue<double>(
+ variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t,
+ x_t.dt);
+
+ porosity_d_data.dphi_dT =
+ dphi_0_dT *
+ (1. + s_therm_exp_data.thermal_volume_strain - biot() * div_u) -
+ (1. - phi_0) * s_therm_exp_data.beta_T_SR;
+}
+
+template struct PorosityModelNonConstantSolidPhaseVolumeFraction<2>;
+template struct PorosityModelNonConstantSolidPhaseVolumeFraction<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h
new file mode 100644
index 0000000000000000000000000000000000000000..f7e0d5d8f09a428005779f300cda2596ac9bd9f1
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h
@@ -0,0 +1,65 @@
+/**
+ * \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
+ */
+
+#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 PorosityDerivativeData
+{
+ double dphi_dT = nan;
+};
+
+struct PorosityData
+{
+ double phi = nan;
+
+ static auto reflect()
+ {
+ using Self = PorosityData;
+ namespace R = ProcessLib::Reflection;
+
+ return std::tuple{R::makeReflectionData("porosity", &Self::phi)};
+ }
+};
+
+struct PorosityModel
+{
+ void eval(SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ PorosityData& porosity_data,
+ PorosityDerivativeData& porosity_d_data) const;
+};
+
+template <int DisplacementDim>
+struct PorosityModelNonConstantSolidPhaseVolumeFraction
+{
+ void eval(
+ SpaceTimeData const& x_t,
+ MediaData const& media_data,
+ BiotData const& biot,
+ StrainData<DisplacementDim> const& strain_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ PorosityData& porosity_data,
+ PorosityDerivativeData& porosity_d_data) const;
+};
+
+extern template struct PorosityModelNonConstantSolidPhaseVolumeFraction<2>;
+extern template struct PorosityModelNonConstantSolidPhaseVolumeFraction<3>;
+
+} // 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
new file mode 100644
index 0000000000000000000000000000000000000000..ad14ab20ff85c3de00210a6740324b59e3866e4f
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h
@@ -0,0 +1,67 @@
+/**
+ * \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
+ */
+
+#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;
+
+ static auto reflect()
+ {
+ using Self = SolidDensityData;
+ namespace R = ProcessLib::Reflection;
+
+ return std::tuple{
+ 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/ConstitutiveRelations/SolidThermalExpansion.cpp b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.cpp
index 39ed38dd266a09329305829b90301a04d6759880..196d01b7904748f5c235d4f271d0399569d44004 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.cpp
@@ -18,6 +18,7 @@ namespace ConstitutiveRelations
template <int DisplacementDim>
void SolidThermalExpansionModel<DisplacementDim>::eval(
SpaceTimeData const& x_t, MediaData const& media_data,
+ TemperatureData const& T_data, ReferenceTemperatureData T0,
SolidThermalExpansionData<DisplacementDim>& out) const
{
namespace MPL = MaterialPropertyLib;
@@ -34,6 +35,9 @@ void SolidThermalExpansionModel<DisplacementDim>::eval(
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
out.beta_T_SR =
Invariants::trace(out.solid_linear_thermal_expansivity_vector);
+
+ double const delta_T = T_data.T - T0();
+ out.thermal_volume_strain = out.beta_T_SR * delta_T;
}
template struct SolidThermalExpansionModel<2>;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h
index 6b7d32057121bd6149b6a83e991ec6aa14cb414d..9f225f30f7d2691173ea2bdde69fd1f136db19af 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidThermalExpansion.h
@@ -21,12 +21,14 @@ struct SolidThermalExpansionData
KelvinVector<DisplacementDim> solid_linear_thermal_expansivity_vector;
double beta_T_SR; /// Isotropic solid phase volumetric thermal expansion
/// coefficient.
+ double thermal_volume_strain = nan;
};
template <int DisplacementDim>
struct SolidThermalExpansionModel
{
void eval(SpaceTimeData const& x_t, MediaData const& media_data,
+ TemperatureData const& T_data, ReferenceTemperatureData T0,
SolidThermalExpansionData<DisplacementDim>& out) const;
};
diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h
index 288d0bbc109bb0eddf2816173c6f94eec410f710..ec9de5fdb91b1fc2a3efbb11e14b5b8828cfd52c 100644
--- a/ProcessLib/TH2M/IntegrationPointData.h
+++ b/ProcessLib/TH2M/IntegrationPointData.h
@@ -36,9 +36,6 @@ struct IntegrationPointData final
typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
- // phase intrinsic densities
- double rhoSR = std::numeric_limits<double>::quiet_NaN();
-
// phase enthalpies
double rho_G_h_G = std::numeric_limits<double>::quiet_NaN();
double rho_L_h_L = std::numeric_limits<double>::quiet_NaN();
@@ -51,10 +48,6 @@ struct IntegrationPointData final
double rho_u_eff = std::numeric_limits<double>::quiet_NaN();
double rho_u_eff_prev = std::numeric_limits<double>::quiet_NaN();
- // porosity
- double phi = std::numeric_limits<double>::quiet_NaN();
- double dphi_dT = std::numeric_limits<double>::quiet_NaN();
-
GlobalDimMatrixType lambda;
GlobalDimVectorType d_CG;
GlobalDimVectorType d_WG;
@@ -64,8 +57,6 @@ struct IntegrationPointData final
GlobalDimVectorType w_GS;
GlobalDimVectorType w_LS;
- double thermal_volume_strain = std::numeric_limits<double>::quiet_NaN();
-
double integration_weight = std::numeric_limits<double>::quiet_NaN();
void pushBackState() { rho_u_eff_prev = rho_u_eff; }
diff --git a/ProcessLib/TH2M/LocalAssemblerInterface.h b/ProcessLib/TH2M/LocalAssemblerInterface.h
index 2efe926f23ebb48a49fc71e8838f58820c7deda4..9e8d86a0863b988bedb2a1a481c85bd38837a952 100644
--- a/ProcessLib/TH2M/LocalAssemblerInterface.h
+++ b/ProcessLib/TH2M/LocalAssemblerInterface.h
@@ -102,18 +102,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
- virtual std::vector<double> const& getIntPtSolidDensity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
- virtual std::vector<double> const& getIntPtPorosity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
virtual std::vector<double> const& getIntPtEnthalpySolid(
const double t,
std::vector<GlobalVector*> const& x,
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 860c34b15f0560bbd2dfe841047a14da1a302e0f..dbd634fb9cc83ebb0ba34e828cc2f127d437585d 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -160,6 +160,8 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ConstitutiveRelations::TemperatureData const T_data{T, T_prev};
ConstitutiveRelations::GasPressureData const pGR_data{pGR};
ConstitutiveRelations::CapillaryPressureData const pCap_data{pCap};
+ ConstitutiveRelations::ReferenceTemperatureData const T0{
+ this->process_data_.reference_temperature(t, pos)[0]};
double const pLR = pGR - pCap;
GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
@@ -196,7 +198,7 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
current_state.swelling_data, ip_cv.swelling_data);
// solid phase linear thermal expansion coefficient
- models.s_therm_exp_model.eval({pos, t, dt}, media_data,
+ models.s_therm_exp_model.eval({pos, t, dt}, media_data, T_data, T0,
ip_cv.s_therm_exp_data);
models.mechanical_strain_model.eval(
@@ -230,6 +232,20 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_out.vapour_pressure_data, current_state.constituent_density_data,
ip_cv.phase_transition_data);
+ models.porosity_model.eval({pos, t, dt}, media_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.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;
@@ -243,48 +259,15 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
vars.liquid_saturation = current_state.S_L_data.S_L;
vars_prev.liquid_saturation = prev_state.S_L_data->S_L;
- auto const rho_ref_SR =
- solid_phase.property(MPL::PropertyType::density)
- .template value<double>(
- vars, pos, t, std::numeric_limits<double>::quiet_NaN());
-
- double const T0 = this->process_data_.reference_temperature(t, pos)[0];
- double const delta_T(T - T0);
- ip_data.thermal_volume_strain =
- ip_cv.s_therm_exp_data.beta_T_SR * delta_T;
-
- // initial porosity
- auto const phi_0 = medium.property(MPL::PropertyType::porosity)
- .template value<double>(vars, pos, t, dt);
-
- auto const phi_S_0 = 1. - phi_0;
-
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- auto const& m = Invariants::identity2;
- double const div_u = m.transpose() * eps;
-
- const double phi_S = phi_S_0 * (1. + ip_data.thermal_volume_strain -
- ip_cv.biot_data() * div_u);
-#else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- const double phi_S = phi_S_0;
-#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
-
- // porosity
- ip_data.phi = 1. - phi_S;
- vars.porosity = ip_data.phi;
-
- // solid phase density
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- auto const rhoSR = rho_ref_SR * (1. - ip_data.thermal_volume_strain +
- (ip_cv.biot_data() - 1.) * div_u);
-#else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- auto const rhoSR = rho_ref_SR;
-#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
+ vars.porosity = ip_out.porosity_data.phi;
auto const& c = ip_cv.phase_transition_data;
- auto const phi_L = current_state.S_L_data.S_L * ip_data.phi;
- auto const phi_G = (1. - current_state.S_L_data.S_L) * ip_data.phi;
+ auto const phi_L =
+ current_state.S_L_data.S_L * ip_out.porosity_data.phi;
+ auto const phi_G =
+ (1. - current_state.S_L_data.S_L) * ip_out.porosity_data.phi;
+ double const phi_S = 1. - ip_out.porosity_data.phi;
// thermal conductivity
ip_data.lambda = MaterialPropertyLib::formEigenTensor<DisplacementDim>(
@@ -301,15 +284,14 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_data.rho_u_eff = phi_G * ip_out.fluid_density_data.rho_GR * c.uG +
phi_L * ip_out.fluid_density_data.rho_LR * c.uL +
- phi_S * rhoSR * u_S;
+ phi_S * ip_out.solid_density_data.rho_SR * u_S;
ip_data.rho_G_h_G =
phi_G * ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G;
ip_data.rho_L_h_L =
phi_L * ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L;
- ip_data.rho_S_h_S = phi_S * rhoSR * ip_data.h_S;
-
- ip_data.rhoSR = rhoSR;
+ ip_data.rho_S_h_S =
+ phi_S * ip_out.solid_density_data.rho_SR * ip_data.h_S;
// for variable output
auto const xmCL = 1. - ip_out.mass_mole_fractions_data.xmWL;
@@ -359,47 +341,27 @@ 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_data.thermal_volume_strain +
- (ip_cv.biot_data() - 1.) * div_u) -
- rho_ref_SR * ip_cv.s_therm_exp_data.beta_T_SR
-#endif
- ;
-
- // porosity
- auto const dphi_0_dT =
- medium[MPL::PropertyType::porosity].template dValue<double>(
- vars, MPL::Variable::temperature, pos, t, dt);
-
- auto const dphi_S_0_dT = -dphi_0_dT;
- const double dphi_S_dT = dphi_S_0_dT
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- * (1. + ip_data.thermal_volume_strain -
- ip_cv.biot_data() * div_u) +
- phi_S_0 * 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 * rhoSR * cpS +
- dphi_S_dT * rhoSR * 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_dT * ip_out.solid_density_data.rho_SR *
+ u_S;
// ds_L_dp_GR = 0;
// ds_G_dp_GR = -ds_L_dp_GR;
double const ds_G_dp_cap = -ip_cv.dS_L_dp_cap();
- // dphi_G_dp_GR = -ds_L_dp_GR * ip_data.phi = 0;
- double const dphi_G_dp_cap = -ip_cv.dS_L_dp_cap() * ip_data.phi;
- // dphi_L_dp_GR = ds_L_dp_GR * ip_data.phi = 0;
- double const dphi_L_dp_cap = ip_cv.dS_L_dp_cap() * ip_data.phi;
+ // dphi_G_dp_GR = -ds_L_dp_GR * ip_out.porosity_data.phi = 0;
+ double const dphi_G_dp_cap =
+ -ip_cv.dS_L_dp_cap() * ip_out.porosity_data.phi;
+ // dphi_L_dp_GR = ds_L_dp_GR * ip_out.porosity_data.phi = 0;
+ double const dphi_L_dp_cap =
+ ip_cv.dS_L_dp_cap() * ip_out.porosity_data.phi;
auto const lambdaGR =
gas_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
@@ -450,7 +412,8 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR;
ip_cv.dlambda_dT = phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT +
- phi_S * dlambda_SR_dT + dphi_S_dT * lambdaSR;
+ phi_S * dlambda_SR_dT -
+ ip_cv.porosity_d_data.dphi_dT * lambdaSR;
// From p_LR = p_GR - p_cap it follows for
// drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR
@@ -486,9 +449,11 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
dphi_L_dT * ip_out.fluid_density_data.rho_LR *
ip_out.enthalpy_data.h_L +
phi_L * drho_LR_dT * ip_out.enthalpy_data.h_L +
- phi_L * ip_out.fluid_density_data.rho_LR * c.dh_L_dT +
- dphi_S_dT * rhoSR * ip_data.h_S + phi_S * drho_SR_dT * ip_data.h_S +
- phi_S * rhoSR * cpS;
+ phi_L * ip_out.fluid_density_data.rho_LR * c.dh_L_dT -
+ ip_cv.porosity_d_data.dphi_dT * ip_out.solid_density_data.rho_SR *
+ 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 =
/*(dphi_G_dp_GR = 0) * ip_out.fluid_density_data.rho_GR * c.uG +*/
@@ -611,33 +576,32 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
/*(ds_L_dp_GR = 0) * current_state.constituent_density_data.rho_C_LR
+*/
s_L * c.drho_C_LR_dp_GR;
- ip_cv.dfC_4_MCpG_dp_GR = drho_C_FR_dp_GR *
- (ip_cv.biot_data() - ip_data.phi) *
- ip_cv.beta_p_SR();
+ ip_cv.dfC_4_MCpG_dp_GR =
+ drho_C_FR_dp_GR * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
+ ip_cv.beta_p_SR();
double const drho_C_FR_dT = s_G * c.drho_C_GR_dT + s_L * c.drho_C_LR_dT;
ip_cv.dfC_4_MCpG_dT =
- drho_C_FR_dT * (ip_cv.biot_data() - ip_data.phi) * ip_cv.beta_p_SR()
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- - rho_C_FR * ip_data.dphi_dT * ip_cv.beta_p_SR()
-#endif
- ;
+ drho_C_FR_dT * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
+ ip_cv.beta_p_SR() -
+ rho_C_FR * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR();
- ip_cv.dfC_4_MCT_dT = drho_C_FR_dT * (ip_cv.biot_data() - ip_data.phi) *
- ip_cv.s_therm_exp_data.beta_T_SR
+ ip_cv.dfC_4_MCT_dT =
+ drho_C_FR_dT * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
+ ip_cv.s_therm_exp_data.beta_T_SR
#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- +
- rho_C_FR * (ip_cv.biot_data() - ip_data.dphi_dT) *
- ip_cv.s_therm_exp_data.beta_T_SR
+ + rho_C_FR * (ip_cv.biot_data() - ip_cv.porosity_d_data.dphi_dT) *
+ ip_cv.s_therm_exp_data.beta_T_SR
#endif
;
ip_cv.dfC_4_MCu_dT = drho_C_FR_dT * ip_cv.biot_data();
- ip_cv.dfC_2a_dp_GR = -ip_data.phi * c.drho_C_GR_dp_GR -
- drho_C_FR_dp_GR * pCap *
- (ip_cv.biot_data() - ip_data.phi) *
- ip_cv.beta_p_SR();
+ ip_cv.dfC_2a_dp_GR =
+ -ip_out.porosity_data.phi * c.drho_C_GR_dp_GR -
+ drho_C_FR_dp_GR * pCap *
+ (ip_cv.biot_data() - ip_out.porosity_data.phi) *
+ ip_cv.beta_p_SR();
double const drho_C_FR_dp_cap =
ds_G_dp_cap * current_state.constituent_density_data.rho_C_GR +
@@ -647,24 +611,22 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
s_L * c.drho_C_LR_dp_LR;
ip_cv.dfC_2a_dp_cap =
- ip_data.phi * (-c.drho_C_LR_dp_LR - drho_C_GR_dp_cap) -
- drho_C_FR_dp_cap * pCap * (ip_cv.biot_data() - ip_data.phi) *
+ ip_out.porosity_data.phi * (-c.drho_C_LR_dp_LR - drho_C_GR_dp_cap) -
+ drho_C_FR_dp_cap * pCap *
+ (ip_cv.biot_data() - ip_out.porosity_data.phi) *
ip_cv.beta_p_SR() +
- rho_C_FR * (ip_cv.biot_data() - ip_data.phi) * ip_cv.beta_p_SR();
+ rho_C_FR * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
+ ip_cv.beta_p_SR();
ip_cv.dfC_2a_dT =
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- ip_data.dphi_dT *
+ ip_cv.porosity_d_data.dphi_dT *
(current_state.constituent_density_data.rho_C_LR -
current_state.constituent_density_data.rho_C_GR) +
-#endif
- ip_data.phi * (c.drho_C_LR_dT - c.drho_C_GR_dT) -
- drho_C_FR_dT * pCap * (ip_cv.biot_data() - ip_data.phi) *
- ip_cv.beta_p_SR()
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- + rho_C_FR * pCap * ip_data.dphi_dT * ip_cv.beta_p_SR()
-#endif
- ;
+ ip_out.porosity_data.phi * (c.drho_C_LR_dT - c.drho_C_GR_dT) -
+ drho_C_FR_dT * pCap *
+ (ip_cv.biot_data() - ip_out.porosity_data.phi) *
+ ip_cv.beta_p_SR() +
+ rho_C_FR * pCap * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR();
ip_cv.dadvection_C_dp_GR = c.drho_C_GR_dp_GR * k_over_mu_G
// + rhoCGR * (dk_over_mu_G_dp_GR = 0)
@@ -698,31 +660,29 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
double const drho_W_FR_dT = s_G * c.drho_W_GR_dT + s_L * c.drho_W_LR_dT;
ip_cv.dfW_2a_dp_GR =
- ip_data.phi * (c.drho_W_LR_dp_GR - c.drho_W_GR_dp_GR);
+ ip_out.porosity_data.phi * (c.drho_W_LR_dp_GR - c.drho_W_GR_dp_GR);
ip_cv.dfW_2b_dp_GR = drho_W_FR_dp_GR * pCap *
- (ip_cv.biot_data() - ip_data.phi) *
+ (ip_cv.biot_data() - ip_out.porosity_data.phi) *
ip_cv.beta_p_SR();
- ip_cv.dfW_2a_dp_cap =
- ip_data.phi * (-c.drho_W_LR_dp_LR - c.drho_W_GR_dp_cap);
+ ip_cv.dfW_2a_dp_cap = ip_out.porosity_data.phi *
+ (-c.drho_W_LR_dp_LR - c.drho_W_GR_dp_cap);
ip_cv.dfW_2b_dp_cap =
- drho_W_FR_dp_cap * pCap * (ip_cv.biot_data() - ip_data.phi) *
+ drho_W_FR_dp_cap * pCap *
+ (ip_cv.biot_data() - ip_out.porosity_data.phi) *
ip_cv.beta_p_SR() +
- rho_W_FR * (ip_cv.biot_data() - ip_data.phi) * ip_cv.beta_p_SR();
+ rho_W_FR * (ip_cv.biot_data() - ip_out.porosity_data.phi) *
+ ip_cv.beta_p_SR();
ip_cv.dfW_2a_dT =
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- ip_data.dphi_dT *
+ ip_cv.porosity_d_data.dphi_dT *
(current_state.rho_W_LR() -
current_state.constituent_density_data.rho_W_GR) +
-#endif
- ip_data.phi * (c.drho_W_LR_dT - c.drho_W_GR_dT);
+ ip_out.porosity_data.phi * (c.drho_W_LR_dT - c.drho_W_GR_dT);
ip_cv.dfW_2b_dT =
- drho_W_FR_dT * pCap * (ip_cv.biot_data() - ip_data.phi) *
- ip_cv.beta_p_SR()
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- - rho_W_FR * pCap * ip_data.dphi_dT * ip_cv.beta_p_SR()
-#endif
- ;
+ drho_W_FR_dT * pCap *
+ (ip_cv.biot_data() - ip_out.porosity_data.phi) *
+ ip_cv.beta_p_SR() -
+ rho_W_FR * pCap * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR();
if (dt == 0.)
{
@@ -1204,22 +1164,17 @@ void TH2MLocalAssembler<
auto const& b = this->process_data_.specific_body_force;
- // porosity
- auto& phi = ip.phi;
-
// volume fraction
- auto const phi_G = s_G * phi;
- auto const phi_L = s_L * phi;
- auto const phi_S = 1. - phi;
-
- // solid phase density
- auto& rho_SR = ip.rhoSR;
+ auto const phi_G = s_G * ip_out.porosity_data.phi;
+ auto const phi_L = s_L * ip_out.porosity_data.phi;
+ auto const phi_S = 1. - ip_out.porosity_data.phi;
auto const rhoGR = ip_out.fluid_density_data.rho_GR;
auto const rhoLR = ip_out.fluid_density_data.rho_LR;
// effective density
- auto const rho = phi_G * rhoGR + phi_L * rhoLR + phi_S * rho_SR;
+ auto const rho = phi_G * rhoGR + phi_L * rhoLR +
+ phi_S * ip_out.solid_density_data.rho_SR;
// abbreviations
const double rho_C_FR =
@@ -1259,9 +1214,12 @@ void TH2MLocalAssembler<
// C-component equation
// ---------------------------------------------------------------------
- MCpG.noalias() += NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * Np * w;
- MCpC.noalias() -=
- NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
+ MCpG.noalias() += NpT * rho_C_FR *
+ (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
+ Np * w;
+ MCpC.noalias() -= NpT * rho_C_FR *
+ (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
+ s_L * Np * w;
if (this->process_data_.apply_mass_lumping)
{
@@ -1269,14 +1227,17 @@ void TH2MLocalAssembler<
{
MCpC.noalias() +=
NpT *
- (phi * (current_state.constituent_density_data.rho_C_LR -
- current_state.constituent_density_data.rho_C_GR) -
- rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
+ (ip_out.porosity_data.phi *
+ (current_state.constituent_density_data.rho_C_LR -
+ current_state.constituent_density_data.rho_C_GR) -
+ rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
+ beta_p_SR) *
s_L_dot * dt / (pCap - pCap_prev) * Np * w;
}
}
- MCT.noalias() -= NpT * rho_C_FR * (alpha_B - phi) * beta_T_SR * Np * w;
+ MCT.noalias() -= NpT * rho_C_FR * (alpha_B - ip_out.porosity_data.phi) *
+ beta_T_SR * Np * w;
MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
using DisplacementDimMatrix =
@@ -1326,22 +1287,27 @@ void TH2MLocalAssembler<
{
fC.noalias() -=
NpT *
- (phi * (current_state.constituent_density_data.rho_C_LR -
- current_state.constituent_density_data.rho_C_GR) -
- rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
+ (ip_out.porosity_data.phi *
+ (current_state.constituent_density_data.rho_C_LR -
+ current_state.constituent_density_data.rho_C_GR) -
+ rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
+ beta_p_SR) *
s_L_dot * w;
}
// fC_III
- fC.noalias() -=
- NpT * phi * (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w;
+ fC.noalias() -= NpT * ip_out.porosity_data.phi *
+ (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w;
// ---------------------------------------------------------------------
// W-component equation
// ---------------------------------------------------------------------
- MWpG.noalias() += NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * Np * w;
- MWpC.noalias() -=
- NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
+ MWpG.noalias() += NpT * rho_W_FR *
+ (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
+ Np * w;
+ MWpC.noalias() -= NpT * rho_W_FR *
+ (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
+ s_L * Np * w;
if (this->process_data_.apply_mass_lumping)
{
@@ -1349,14 +1315,17 @@ void TH2MLocalAssembler<
{
MWpC.noalias() +=
NpT *
- (phi * (current_state.rho_W_LR() -
- current_state.constituent_density_data.rho_W_GR) -
- rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
+ (ip_out.porosity_data.phi *
+ (current_state.rho_W_LR() -
+ current_state.constituent_density_data.rho_W_GR) -
+ rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
+ beta_p_SR) *
s_L_dot * dt / (pCap - pCap_prev) * Np * w;
}
}
- MWT.noalias() -= NpT * rho_W_FR * (alpha_B - phi) * beta_T_SR * Np * w;
+ MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
+ beta_T_SR * Np * w;
MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
@@ -1404,14 +1373,16 @@ void TH2MLocalAssembler<
{
fW.noalias() -=
NpT *
- (phi * (current_state.rho_W_LR() -
- current_state.constituent_density_data.rho_W_GR) -
- rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
+ (ip_out.porosity_data.phi *
+ (current_state.rho_W_LR() -
+ current_state.constituent_density_data.rho_W_GR) -
+ rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
+ beta_p_SR) *
s_L_dot * w;
}
- fW.noalias() -=
- NpT * phi * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
+ fW.noalias() -= NpT * ip_out.porosity_data.phi *
+ (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
// ---------------------------------------------------------------------
// - temperature equation
@@ -1684,22 +1655,17 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
auto const& b = this->process_data_.specific_body_force;
- // porosity
- auto& phi = ip.phi;
-
// volume fraction
- auto const phi_G = s_G * phi;
- auto const phi_L = s_L * phi;
- auto const phi_S = 1. - phi;
-
- // solid phase density
- auto& rho_SR = ip.rhoSR;
+ auto const phi_G = s_G * ip_out.porosity_data.phi;
+ auto const phi_L = s_L * ip_out.porosity_data.phi;
+ auto const phi_S = 1. - ip_out.porosity_data.phi;
auto const rhoGR = ip_out.fluid_density_data.rho_GR;
auto const rhoLR = ip_out.fluid_density_data.rho_LR;
// effective density
- auto const rho = phi_G * rhoGR + phi_L * rhoLR + phi_S * rho_SR;
+ auto const rho = phi_G * rhoGR + phi_L * rhoLR +
+ phi_S * ip_out.solid_density_data.rho_SR;
// abbreviations
const double rho_C_FR =
@@ -1739,9 +1705,12 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
// C-component equation
// ---------------------------------------------------------------------
- MCpG.noalias() += NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * Np * w;
- MCpC.noalias() -=
- NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
+ MCpG.noalias() += NpT * rho_C_FR *
+ (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
+ Np * w;
+ MCpC.noalias() -= NpT * rho_C_FR *
+ (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
+ s_L * Np * w;
if (this->process_data_.apply_mass_lumping)
{
@@ -1749,14 +1718,17 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
{
MCpC.noalias() +=
NpT *
- (phi * (current_state.constituent_density_data.rho_C_LR -
- current_state.constituent_density_data.rho_C_GR) -
- rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
+ (ip_out.porosity_data.phi *
+ (current_state.constituent_density_data.rho_C_LR -
+ current_state.constituent_density_data.rho_C_GR) -
+ rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
+ beta_p_SR) *
s_L_dot * dt / (pCap - pCap_prev) * Np * w;
}
}
- MCT.noalias() -= NpT * rho_C_FR * (alpha_B - phi) * beta_T_SR * Np * w;
+ MCT.noalias() -= NpT * rho_C_FR * (alpha_B - ip_out.porosity_data.phi) *
+ beta_T_SR * Np * w;
// d (fC_4_MCT * T_dot)/d T
local_Jac
.template block<C_size, temperature_size>(C_index,
@@ -1863,9 +1835,11 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
{
// fC_2 = \int a * s_L_dot
auto const a =
- phi * (current_state.constituent_density_data.rho_C_LR -
- current_state.constituent_density_data.rho_C_GR) -
- rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR;
+ ip_out.porosity_data.phi *
+ (current_state.constituent_density_data.rho_C_LR -
+ current_state.constituent_density_data.rho_C_GR) -
+ rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
+ beta_p_SR;
fC.noalias() -= NpT * a * s_L_dot * w;
local_Jac.template block<C_size, C_size>(C_index, C_index)
@@ -1888,32 +1862,34 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
{
// fC_3 = \int phi * a
double const a = s_G * rho_C_GR_dot + s_L * rho_C_LR_dot;
- fC.noalias() -= NpT * phi * a * w;
+ fC.noalias() -= NpT * ip_out.porosity_data.phi * a * w;
local_Jac.template block<C_size, C_size>(C_index, C_index)
- .noalias() += NpT * phi * ip_cv.dfC_3a_dp_GR * Np * w;
+ .noalias() +=
+ NpT * ip_out.porosity_data.phi * ip_cv.dfC_3a_dp_GR * Np * w;
local_Jac.template block<C_size, W_size>(C_index, W_index)
- .noalias() += NpT * phi * ip_cv.dfC_3a_dp_cap * Np * w;
+ .noalias() +=
+ NpT * ip_out.porosity_data.phi * ip_cv.dfC_3a_dp_cap * Np * w;
local_Jac
.template block<C_size, temperature_size>(C_index,
temperature_index)
.noalias() += NpT *
- (
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- ip.dphi_dT * a +
-#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- phi * ip_cv.dfC_3a_dT) *
+ (ip_cv.porosity_d_data.dphi_dT * a +
+ ip_out.porosity_data.phi * ip_cv.dfC_3a_dT) *
NT * w;
}
// ---------------------------------------------------------------------
// W-component equation
// ---------------------------------------------------------------------
- MWpG.noalias() += NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * Np * w;
- MWpC.noalias() -=
- NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
+ MWpG.noalias() += NpT * rho_W_FR *
+ (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
+ Np * w;
+ MWpC.noalias() -= NpT * rho_W_FR *
+ (alpha_B - ip_out.porosity_data.phi) * beta_p_SR *
+ s_L * Np * w;
if (this->process_data_.apply_mass_lumping)
{
@@ -1921,14 +1897,17 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
{
MWpC.noalias() +=
NpT *
- (phi * (current_state.rho_W_LR() -
- current_state.constituent_density_data.rho_W_GR) -
- rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
+ (ip_out.porosity_data.phi *
+ (current_state.rho_W_LR() -
+ current_state.constituent_density_data.rho_W_GR) -
+ rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) *
+ beta_p_SR) *
s_L_dot * dt / (pCap - pCap_prev) * Np * w;
}
}
- MWT.noalias() -= NpT * rho_W_FR * (alpha_B - phi) * beta_T_SR * Np * w;
+ MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
+ beta_T_SR * Np * w;
MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
@@ -1997,10 +1976,11 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
// fW_2 = \int (f - g) * s_L_dot
if (!this->process_data_.apply_mass_lumping)
{
- double const f =
- phi * (current_state.rho_W_LR() -
- current_state.constituent_density_data.rho_W_GR);
- double const g = rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR;
+ double const f = ip_out.porosity_data.phi *
+ (current_state.rho_W_LR() -
+ current_state.constituent_density_data.rho_W_GR);
+ double const g = rho_W_FR * pCap *
+ (alpha_B - ip_out.porosity_data.phi) * beta_p_SR;
fW.noalias() -= NpT * (f - g) * s_L_dot * w;
@@ -2026,26 +2006,23 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
}
// fW_3 = \int phi * a
- fW.noalias() -=
- NpT * phi * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
+ fW.noalias() -= NpT * ip_out.porosity_data.phi *
+ (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
- NpT * phi * ip_cv.dfW_3a_dp_GR * Np * w;
+ NpT * ip_out.porosity_data.phi * ip_cv.dfW_3a_dp_GR * Np * w;
local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
- NpT * phi * ip_cv.dfW_3a_dp_cap * Np * w;
+ NpT * ip_out.porosity_data.phi * ip_cv.dfW_3a_dp_cap * Np * w;
local_Jac
.template block<W_size, temperature_size>(W_index,
temperature_index)
- .noalias() +=
- NpT *
- (
-#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- ip.dphi_dT * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) +
-#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- phi * ip_cv.dfW_3a_dT) *
- NT * w;
+ .noalias() += NpT *
+ (ip_cv.porosity_d_data.dphi_dT *
+ (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) +
+ ip_out.porosity_data.phi * ip_cv.dfW_3a_dT) *
+ NT * w;
// ---------------------------------------------------------------------
// - temperature equation
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index 3a5a1a4e25f7570267f61023e2c06e199a5a4692..da9643ed1706c890f62afe42056776b839226845 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -235,26 +235,6 @@ private:
// be stored in some container to avoid defining one method for each
// variable.
- virtual std::vector<double> const& getIntPtSolidDensity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const override
- {
- return ProcessLib::getIntegrationPointScalarData(_ip_data,
- &IpData::rhoSR, cache);
- }
-
- virtual std::vector<double> const& getIntPtPorosity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const override
- {
- return ProcessLib::getIntegrationPointScalarData(_ip_data, &IpData::phi,
- cache);
- }
-
virtual std::vector<double> const& getIntPtEnthalpySolid(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index b5691b3012a74f9eb676b0c11bc0e9c439eb716e..34d2e0d1db4de1e5c312f543bbb073b56b3420b2 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -186,13 +186,6 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
&LocalAssemblerInterface<
DisplacementDim>::getIntPtDiffusionVelocityLiquidLiquid);
- add_secondary_variable(
- "porosity", 1,
- &LocalAssemblerInterface<DisplacementDim>::getIntPtPorosity);
- add_secondary_variable(
- "solid_density", 1,
- &LocalAssemblerInterface<DisplacementDim>::getIntPtSolidDensity);
-
add_secondary_variable(
"enthalpy_solid", 1,
&LocalAssemblerInterface<DisplacementDim>::getIntPtEnthalpySolid);
diff --git a/Tests/Data/TH2M/H2/dissolution_diffusion/bourgeat.prj b/Tests/Data/TH2M/H2/dissolution_diffusion/bourgeat.prj
index c3a73319ad2edf4961c49ff9aa2474844c37ad07..f07c0c3fd0dc5f22882eee6c2288feec927e0c18 100644
--- a/Tests/Data/TH2M/H2/dissolution_diffusion/bourgeat.prj
+++ b/Tests/Data/TH2M/H2/dissolution_diffusion/bourgeat.prj
@@ -522,7 +522,7 @@
<vtkdiff>
<regex>result_bourgeat_ts_.*.vtu</regex>
<field>gas_pressure</field>
- <absolute_tolerance>5.6e-3</absolute_tolerance>
+ <absolute_tolerance>8e-3</absolute_tolerance>
<relative_tolerance>0</relative_tolerance>
</vtkdiff>
<vtkdiff>
@@ -576,7 +576,7 @@
<vtkdiff>
<regex>result_bourgeat_ts_.*.vtu</regex>
<field>liquid_density</field>
- <absolute_tolerance>1e-10</absolute_tolerance>
+ <absolute_tolerance>1.3e-10</absolute_tolerance>
<relative_tolerance>0</relative_tolerance>
</vtkdiff>
<vtkdiff>
diff --git a/Tests/Data/TH2M/H2/dissolution_diffusion/continuous_injection.prj b/Tests/Data/TH2M/H2/dissolution_diffusion/continuous_injection.prj
index 2ed7a815b54f4018b931af3ac97b1cf11524b800..b5905dda6e01fd6f882d970fda533e1a1c801bb5 100644
--- a/Tests/Data/TH2M/H2/dissolution_diffusion/continuous_injection.prj
+++ b/Tests/Data/TH2M/H2/dissolution_diffusion/continuous_injection.prj
@@ -259,7 +259,7 @@
<property>
<name>porosity</name>
<type>Constant</type>
- <value>0.15</value>
+ <value>0.15000000000000002</value>
</property>
<property>
<name>tortuosity</name>
diff --git a/Tests/Data/TH2M/TH/Ogata-Banks/ogata-banks.prj b/Tests/Data/TH2M/TH/Ogata-Banks/ogata-banks.prj
index 6ae4bacbc620b4fe2f53925358a9a8f222d704a4..724c5ed5e0ada6702d3eb543bb6de09a0eccbd1a 100644
--- a/Tests/Data/TH2M/TH/Ogata-Banks/ogata-banks.prj
+++ b/Tests/Data/TH2M/TH/Ogata-Banks/ogata-banks.prj
@@ -456,8 +456,8 @@
<vtkdiff>
<regex>result_ogata-banks_ts_.*.vtu</regex>
<field>temperature_interpolated</field>
- <absolute_tolerance>1e-09</absolute_tolerance>
+ <absolute_tolerance>1.7e-09</absolute_tolerance>
<relative_tolerance>0</relative_tolerance>
</vtkdiff>
</test_definition>
-</OpenGeoSysProject>
\ No newline at end of file
+</OpenGeoSysProject>