diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index 7bcdc10f1174067792fdc668e233eea8b5582312..2cd3860282d84a67feef5e0595dbba2076823ffa 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -41,6 +41,7 @@
#include "TEquation.h"
#include "ThermalConductivity.h"
#include "TotalStress.h"
+#include "UEquation.h"
#include "VapourPartialPressure.h"
#include "Viscosity.h"
#include "WEquation.h"
@@ -186,6 +187,8 @@ struct ConstitutiveTempData
FT1Data fT_1;
FT2Data<DisplacementDim> fT_2;
FT3Data<DisplacementDim> fT_3;
+
+ FU2KUpCData fu_2_KupC;
};
/// Data that stores intermediate values (derivatives), which are not needed
@@ -213,6 +216,7 @@ struct DerivativesData
FW4LWpCDerivativeData<DisplacementDim> dfW_4_LWpC;
FT1DerivativeData dfT_1;
FT2DerivativeData<DisplacementDim> dfT_2;
+ FU2KUpCDerivativeData dfu_2_KupC;
};
} // namespace ConstitutiveRelations
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index 73d5a7c922dd9df0c32f16839f933313b726ef1e..c79d9e81fa825c3ba5890b34f674438d6fbbe7c4 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -33,6 +33,7 @@
#include "TEquation.h"
#include "ThermalConductivity.h"
#include "TotalStress.h"
+#include "UEquation.h"
#include "Viscosity.h"
#include "WEquation.h"
@@ -114,6 +115,8 @@ struct ConstitutiveModels
FT1Model fT_1_model;
FT2Model<DisplacementDim> fT_2_model;
FT3Model<DisplacementDim> fT_3_model;
+
+ FU2KUpCModel fu_2_KupC_model;
};
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..07f391f74d88d7c1534a714d91756e9df518658a
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp
@@ -0,0 +1,34 @@
+/**
+ * \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 "UEquation.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+void FU2KUpCModel::eval(BiotData const& biot_data,
+ BishopsData const& chi_S_L,
+ FU2KUpCData& fu_2_KupC) const
+{
+ fu_2_KupC.m = biot_data() * chi_S_L.chi_S_L;
+}
+
+void FU2KUpCModel::dEval(BiotData const& biot_data,
+ BishopsData const& chi_S_L,
+ CapillaryPressureData const& p_cap,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ FU2KUpCDerivativeData& dfu_2_KupC) const
+{
+ dfu_2_KupC.dp_cap =
+ biot_data() * chi_S_L.dchi_dS_L * dS_L_dp_cap() * p_cap();
+}
+
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h
new file mode 100644
index 0000000000000000000000000000000000000000..bffa79f508b57944cf18655d9a9846c29e57dbbc
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h
@@ -0,0 +1,44 @@
+/**
+ * \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 "Bishops.h"
+#include "Saturation.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+struct FU2KUpCData
+{
+ double m = nan;
+};
+
+struct FU2KUpCDerivativeData
+{
+ double dp_cap = nan;
+};
+
+struct FU2KUpCModel
+{
+ void eval(BiotData const& biot_data,
+ BishopsData const& chi_S_L,
+ FU2KUpCData& fu_2_KupC) const;
+
+ void dEval(BiotData const& biot_data,
+ BishopsData const& chi_S_L,
+ CapillaryPressureData const& p_cap,
+ SaturationDataDeriv const& dS_L_dp_cap,
+ FU2KUpCDerivativeData& dfu_2_KupC) const;
+};
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 31e8e89f01178919c22ec16a0f2de9b250d2d990..2f2fe779e7af872abc27b6367bddc02da4e57cbc 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -473,6 +473,9 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ConstitutiveRelations::SpecificBodyForceData<DisplacementDim>{
this->process_data_.specific_body_force},
ip_cv.fT_3);
+
+ models.fu_2_KupC_model.eval(ip_cv.biot_data, ip_cv.chi_S_L,
+ ip_cv.fu_2_KupC);
}
return {ip_constitutive_data, ip_constitutive_variables};
@@ -540,8 +543,8 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
double const T = NT.dot(temperature);
double const T_prev = NT.dot(temperature_prev);
ConstitutiveRelations::TemperatureData const T_data{T, T_prev};
- double const pCap = Np.dot(capillary_pressure);
- ConstitutiveRelations::CapillaryPressureData const pCap_data{pCap};
+ ConstitutiveRelations::CapillaryPressureData const pCap_data{
+ Np.dot(capillary_pressure)};
models.advection_model.dEval(current_state.constituent_density_data,
ip_out.permeability_data,
@@ -702,6 +705,12 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
this->process_data_.specific_body_force},
ip_cv.viscosity_data,
ip_dd.dfT_2);
+
+ models.fu_2_KupC_model.dEval(ip_cv.biot_data,
+ ip_cv.chi_S_L,
+ pCap_data,
+ ip_cv.dS_L_dp_cap,
+ ip_dd.dfu_2_KupC);
}
return ip_d_data;
@@ -1165,7 +1174,8 @@ void TH2MLocalAssembler<
KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
- KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
+ KUpC.noalias() +=
+ BuT * Invariants::identity2 * Np * (ip_cv.fu_2_KupC.m * w);
fU.noalias() -=
(BuT * current_state.eff_stress_data.sigma -
@@ -1742,15 +1752,16 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
// dfU_2/dp_GR = dKUpG/dp_GR * p_GR + KUpG. The former is zero, the
// latter is handled below.
- KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
+ KUpC.noalias() +=
+ BuT * Invariants::identity2 * Np * (ip_cv.fu_2_KupC.m * w);
// dfU_2/dp_cap = dKUpC/dp_cap * p_cap + KUpC. The former is handled
// here, the latter below.
local_Jac
.template block<displacement_size, W_size>(displacement_index,
W_index)
- .noalias() += BuT * alpha_B * ip_cv.chi_S_L.dchi_dS_L *
- ip_cv.dS_L_dp_cap() * pCap * m * Np * w;
+ .noalias() +=
+ BuT * Invariants::identity2 * Np * (ip_dd.dfu_2_KupC.dp_cap * w);
local_Jac
.template block<displacement_size, displacement_size>(