From 29b62eae9121aa026230f84d7653da366faa6e5d Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Tue, 28 Sep 2021 16:13:07 +0200
Subject: [PATCH] [TH2M] Enable to use pressure, strain, stress dependent
permeability models
---
ProcessLib/TH2M/TH2MFEM-impl.h | 43 +++++++++++++++++++++++++++-------
ProcessLib/TH2M/TH2MFEM.h | 2 ++
2 files changed, 36 insertions(+), 9 deletions(-)
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index d7763bb9713..115b76ed869 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -182,7 +182,41 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
medium.property(MPL::PropertyType::permeability)
.value(vars, pos, t, dt));
+ auto const Bu =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunctionDisplacement::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ gradNu, Nu, x_coord, _is_axially_symmetric);
+
+ auto& eps = ip_data.eps;
+ eps.noalias() = Bu * displacement;
+
// relative permeability
+ // Set mechanical variables for the intrinsic permeability model
+ // For stress dependent permeability.
+ {
+ // Note: if Bishop model is available, ip_data.s_L in the following
+ // computation should be replaced with the Bishop value.
+ auto const sigma_total =
+ (_ip_data[ip].sigma_eff - ip_data.alpha_B *
+ (pGR - ip_data.s_L * pCap) *
+ Invariants::identity2)
+ .eval();
+
+ vars[static_cast<int>(MPL::Variable::total_stress)]
+ .emplace<SymmetricTensor>(
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
+ sigma_total));
+ }
+ // For strain dependent permeability
+ vars[static_cast<int>(MPL::Variable::volumetric_strain)] =
+ Invariants::trace(eps);
+ vars[static_cast<int>(MPL::Variable::equivalent_plastic_strain)] =
+ _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
+ vars[static_cast<int>(MPL::Variable::mechanical_strain)]
+ .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
+ eps);
+
ip_data.k_rel_G =
medium
.property(
@@ -211,15 +245,6 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
dthermal_strain = ip_data.alpha_T_SR * T_dot * dt;
- auto const Bu =
- LinearBMatrix::computeBMatrix<DisplacementDim,
- ShapeFunctionDisplacement::NPOINTS,
- typename BMatricesType::BMatrixType>(
- gradNu, Nu, x_coord, _is_axially_symmetric);
-
- auto& eps = ip_data.eps;
- eps.noalias() = Bu * displacement;
-
auto& eps_prev = ip_data.eps_prev;
auto& eps_m = ip_data.eps_m;
auto& eps_m_prev = ip_data.eps_m_prev;
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index b770cb1b838..a62a3267368 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -70,6 +70,8 @@ public:
static int const KelvinVectorSize =
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+ using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
+
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
TH2MLocalAssembler(TH2MLocalAssembler const&) = delete;
--
GitLab