From e9456935c7dfa86f72af183df8538f945055176f Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon, 8 Jun 2020 16:31:39 +0200
Subject: [PATCH] [RichardsMechanics] use stress dependent permeability
---
.../RichardsMechanicsFEM-impl.h | 45 +++++++++++++++----
.../RichardsMechanics/RichardsMechanicsFEM.h | 2 +
2 files changed, 38 insertions(+), 9 deletions(-)
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index f3549417cde..10eed56d986 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -458,8 +458,18 @@ void RichardsMechanicsLocalAssembler<
double const k_rel =
medium->property(MPL::PropertyType::relative_permeability)
.template value<double>(variables, x_position, t, dt);
- auto const mu = liquid_phase.property(MPL::PropertyType::viscosity)
- .template value<double>(variables, x_position, t, dt);
+ auto const mu =
+ liquid_phase.property(MPL::PropertyType::viscosity)
+ .template value<double>(variables, x_position, t, dt);
+
+ // For stress dependent permeability.
+ variables[static_cast<int>(MPL::Variable::stress)]
+ .emplace<SymmetricTensor>(
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
+ (_ip_data[ip].sigma_eff +
+ alpha * chi_S_L * identity2 * p_cap_ip)
+ .eval()));
+
auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
solid_phase.property(MPL::PropertyType::permeability)
.value(variables, x_position, t, dt));
@@ -510,8 +520,7 @@ void RichardsMechanicsLocalAssembler<
//
K.template block<displacement_size, pressure_size>(displacement_index,
pressure_index)
- .noalias() -=
- B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
+ .noalias() -= B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
//
// pressure equation, displacement part.
@@ -779,8 +788,17 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
double const k_rel =
medium->property(MPL::PropertyType::relative_permeability)
.template value<double>(variables, x_position, t, dt);
- auto const mu = liquid_phase.property(MPL::PropertyType::viscosity)
- .template value<double>(variables, x_position, t, dt);
+ auto const mu =
+ liquid_phase.property(MPL::PropertyType::viscosity)
+ .template value<double>(variables, x_position, t, dt);
+
+ // For stress dependent permeability.
+ variables[static_cast<int>(MPL::Variable::stress)]
+ .emplace<SymmetricTensor>(
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
+ (_ip_data[ip].sigma_eff +
+ alpha * chi_S_L * identity2 * p_cap_ip)
+ .eval()));
auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
solid_phase.property(MPL::PropertyType::permeability)
.value(variables, x_position, t, dt));
@@ -1408,7 +1426,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
typename BMatricesType::BMatrixType>(
dNdx_u, N_u, x_coord, _is_axially_symmetric);
-
double p_cap_ip;
NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
@@ -1527,14 +1544,24 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
porosity;
}
}
- auto const mu = liquid_phase.property(MPL::PropertyType::viscosity)
- .template value<double>(variables, x_position, t, dt);
+ auto const mu =
+ liquid_phase.property(MPL::PropertyType::viscosity)
+ .template value<double>(variables, x_position, t, dt);
auto const rho_LR =
liquid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
+
+ // For stress dependent permeability.
+ variables[static_cast<int>(MPL::Variable::stress)]
+ .emplace<SymmetricTensor>(
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
+ (_ip_data[ip].sigma_eff +
+ alpha * chi_S_L * identity2 * p_cap_ip)
+ .eval()));
auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
solid_phase.property(MPL::PropertyType::permeability)
.value(variables, x_position, t, dt));
+
double const k_rel =
medium->property(MPL::PropertyType::relative_permeability)
.template value<double>(variables, x_position, t, dt);
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
index c4ce33132c0..67fb63b4f85 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -69,6 +69,8 @@ public:
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
+ using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
+
RichardsMechanicsLocalAssembler(RichardsMechanicsLocalAssembler const&) =
delete;
RichardsMechanicsLocalAssembler(RichardsMechanicsLocalAssembler&&) = delete;
--
GitLab