From 0ff7aa0cd60e72843dc7d0cf5c8e1a675a259e96 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Mon, 15 Apr 2024 14:19:43 +0200
Subject: [PATCH] [PL/TH2M] Extract MWT part of W4 equation
---
.../ConstitutiveRelations/ConstitutiveData.h | 1 +
.../ConstitutiveModels.h | 1 +
.../TH2M/ConstitutiveRelations/WEquation.cpp | 22 +++++++++++++++++++
.../TH2M/ConstitutiveRelations/WEquation.h | 21 ++++++++++++++++++
ProcessLib/TH2M/TH2MFEM-impl.h | 17 ++++++++------
5 files changed, 55 insertions(+), 7 deletions(-)
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index 7312f86bc8c..582d8f8215b 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -175,6 +175,7 @@ struct ConstitutiveTempData
FW4LWTData<DisplacementDim> fW_4_LWT;
FW4MWpGData fW_4_MWpG;
FW4MWpCData fW_4_MWpC;
+ FW4MWTData fW_4_MWT;
using DisplacementDimVector = Eigen::Matrix<double, DisplacementDim, 1>;
using DisplacementDimMatrix =
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index 4c4eec3ad79..01b100cbe32 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -101,6 +101,7 @@ struct ConstitutiveModels
FW4LWTModel<DisplacementDim> fW_4_LWT_model;
FW4MWpGModel fW_4_MWpG_model;
FW4MWpCModel fW_4_MWpC_model;
+ FW4MWTModel<DisplacementDim> fW_4_MWT_model;
};
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
index c269bad6fd8..721c36ad78d 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
+++ b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp
@@ -416,5 +416,27 @@ void FW4MWpCModel::eval(BiotData const& biot_data,
(S_L - S_L_data_prev->S_L);
}
+template <int DisplacementDim>
+void FW4MWTModel<DisplacementDim>::eval(
+ BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ FW4MWTData& fW_4_MWT) const
+{
+ auto const S_L = S_L_data.S_L;
+ auto const S_G = 1. - S_L;
+ double const rho_W_FR =
+ S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR();
+
+ fW_4_MWT.m = -rho_W_FR * (biot_data() - porosity_data.phi) *
+ s_therm_exp_data.beta_T_SR;
+}
+
+template struct FW4MWTModel<2>;
+template struct FW4MWTModel<3>;
+
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
index ac2644843b0..e2f99341be5 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h
@@ -239,5 +239,26 @@ struct FW4MWpCModel
FW4MWpCData& fW_4_MWpC) const;
};
+struct FW4MWTData
+{
+ double m = nan;
+};
+
+template <int DisplacementDim>
+struct FW4MWTModel
+{
+ void eval(
+ BiotData const& biot_data,
+ ConstituentDensityData const& constituent_density_data,
+ PorosityData const& porosity_data,
+ PureLiquidDensityData const& rho_W_LR,
+ SaturationData const& S_L_data,
+ SolidThermalExpansionData<DisplacementDim> const& s_therm_exp_data,
+ FW4MWTData& fW_4_MWT) const;
+};
+
+extern template struct FW4MWTModel<2>;
+extern template struct FW4MWTModel<3>;
+
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index c4dc4dc316c..df7288e7eeb 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -411,6 +411,14 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_cv.beta_p_SR,
ip_cv.fW_4_MWpC);
+ models.fW_4_MWT_model.eval(ip_cv.biot_data,
+ current_state.constituent_density_data,
+ ip_out.porosity_data,
+ current_state.rho_W_LR,
+ current_state.S_L_data,
+ ip_cv.s_therm_exp_data,
+ ip_cv.fW_4_MWT);
+
// for variable output
auto const xmCL = 1. - ip_out.mass_mole_fractions_data.xmWL;
@@ -1089,8 +1097,6 @@ void TH2MLocalAssembler<
double const pCap = Np.dot(capillary_pressure);
double const pCap_prev = Np.dot(capillary_pressure_prev);
- double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
-
auto const s_L = current_state.S_L_data.S_L;
auto const s_G = 1. - s_L;
auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
@@ -1166,8 +1172,7 @@ void TH2MLocalAssembler<
}
}
- MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
- beta_T_SR * Np * w;
+ MWT.noalias() += NpT * ip_cv.fW_4_MWT.m * Np * w;
MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
@@ -1451,7 +1456,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
double const pGR_prev = Np.dot(gas_pressure_prev);
double const pCap_prev = Np.dot(capillary_pressure_prev);
double const T_prev = NT.dot(temperature_prev);
- double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR;
auto& s_L = current_state.S_L_data.S_L;
auto const s_G = 1. - s_L;
@@ -1614,8 +1618,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
}
}
- MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) *
- beta_T_SR * Np * w;
+ MWT.noalias() += NpT * ip_cv.fW_4_MWT.m * Np * w;
MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
--
GitLab