diff --git a/Documentation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/t_intrinsic_permeability.md b/Documentation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/t_intrinsic_permeability.md
deleted file mode 100644
index f945fe4d68fe60bacc884d5d1504e0ad52e4f6c7..0000000000000000000000000000000000000000
--- a/Documentation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/t_intrinsic_permeability.md
+++ /dev/null
@@ -1 +0,0 @@
-\copydoc ProcessLib::RichardsMechanics::RichardsMechanicsProcessData::intrinsic_permeability
diff --git a/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
index 40fe159072d5b2284b36c0fc66f4add87757dbb5..44a939c5534a60e2ca3f3d7ff0a4e8be6174b451 100644
--- a/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
@@ -105,15 +105,6 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
parameters, config);
- // Intrinsic permeability
- auto& intrinsic_permeability = findParameter<double>(
- config,
- //! \ogs_file_param_special{prj__processes__process__RICHARDS_MECHANICS__intrinsic_permeability}
- "intrinsic_permeability", parameters, 1);
-
- DBUG("Use '%s' as intrinsic conductivity parameter.",
- intrinsic_permeability.name.c_str());
-
// Fluid bulk modulus
auto& fluid_bulk_modulus = findParameter<double>(
config,
@@ -197,7 +188,6 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
material_ids,
std::move(flow_material),
std::move(solid_constitutive_relations),
- intrinsic_permeability,
fluid_bulk_modulus,
biot_coefficient,
solid_density,
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 250df5ed328f776953e2d48ce266547534bf471a..dc0f58a397f7e664a4285c1fafeb6da13c9b14a1 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -13,6 +13,8 @@
#include "RichardsMechanicsFEM.h"
+#include <cassert>
+
#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
#include "MathLib/KelvinVector.h"
#include "NumLib/Function/Interpolation.h"
@@ -22,6 +24,29 @@ namespace ProcessLib
{
namespace RichardsMechanics
{
+template <int DisplacementDim>
+Eigen::Matrix<double, DisplacementDim, DisplacementDim> intrinsicPermeability(
+ double const t, SpatialPosition const& x_position, int const material_id,
+ RichardsFlow::RichardsFlowMaterialProperties const& material)
+{
+ const Eigen::MatrixXd& permeability =
+ material.getPermeability(material_id, t, x_position, DisplacementDim);
+ if (permeability.rows() == DisplacementDim)
+ {
+ return permeability;
+ }
+ if (permeability.rows() == 1)
+ {
+ return Eigen::Matrix<double, DisplacementDim,
+ DisplacementDim>::Identity() *
+ permeability(0, 0);
+ }
+
+ OGS_FATAL(
+ "Intrinsic permeability dimension is neither %d nor one, but %dx%d.",
+ DisplacementDim, permeability.rows(), permeability.cols());
+}
+
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
@@ -191,12 +216,14 @@ void RichardsMechanicsLocalAssembler<
double const k_rel =
_process_data.flow_material->getRelativePermeability(
t, x_position, -p_cap_ip, temperature, S_L);
- auto const mu = _process_data.flow_material->getFluidViscosity(
+
+ double const mu = _process_data.flow_material->getFluidViscosity(
-p_cap_ip, temperature);
- double const K_intrinsic =
- _process_data.intrinsic_permeability(t, x_position)[0];
- double const K_over_mu = K_intrinsic / mu;
+ GlobalDimMatrixType const rho_K_over_mu =
+ intrinsicPermeability<DisplacementDim>(
+ t, x_position, material_id, *_process_data.flow_material) *
+ (rho_LR * k_rel / mu);
//
// displacement equation, displacement part
@@ -231,11 +258,10 @@ void RichardsMechanicsLocalAssembler<
K.template block<pressure_size, pressure_size>(pressure_index,
pressure_index)
- .noalias() +=
- dNdx_p.transpose() * rho_LR * k_rel * K_over_mu * dNdx_p * w;
+ .noalias() += dNdx_p.transpose() * rho_K_over_mu * dNdx_p * w;
rhs.template segment<pressure_size>(pressure_index).noalias() +=
- dNdx_p.transpose() * rho_LR * rho_LR * k_rel * K_over_mu * b * w;
+ dNdx_p.transpose() * rho_LR * rho_K_over_mu * b * w;
//
// displacement equation, pressure part
@@ -401,12 +427,15 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
double const k_rel =
_process_data.flow_material->getRelativePermeability(
t, x_position, -p_cap_ip, temperature, S_L);
- auto const mu = _process_data.flow_material->getFluidViscosity(
+
+ double const mu = _process_data.flow_material->getFluidViscosity(
-p_cap_ip, temperature);
- double const K_intrinsic =
- _process_data.intrinsic_permeability(t, x_position)[0];
- double const K_over_mu = K_intrinsic / mu;
+ GlobalDimMatrixType const rho_Ki_over_mu =
+ intrinsicPermeability<DisplacementDim>(
+ t, x_position, material_id, *_process_data.flow_material) *
+ (rho_LR / mu);
+
//
// displacement equation, displacement part
//
@@ -461,7 +490,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
// pressure equation, pressure part.
//
laplace_p.noalias() +=
- dNdx_p.transpose() * rho_LR * k_rel * K_over_mu * dNdx_p * w;
+ dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
double const a0 = (alpha - porosity) / K_SR;
double const specific_storage =
@@ -508,17 +537,17 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
local_Jac
.template block<pressure_size, pressure_size>(pressure_index,
pressure_index)
- .noalias() += dNdx_p.transpose() * rho_LR * K_over_mu * grad_p_cap *
+ .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
dk_rel_dS_l * dS_L_dp_cap * N_p * w;
local_Jac
.template block<pressure_size, pressure_size>(pressure_index,
pressure_index)
- .noalias() += dNdx_p.transpose() * rho_LR * rho_LR * K_over_mu * b *
+ .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
dk_rel_dS_l * dS_L_dp_cap * N_p * w;
local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
- dNdx_p.transpose() * rho_LR * rho_LR * k_rel * K_over_mu * b * w;
+ dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
}
// pressure equation, pressure part.
@@ -626,6 +655,9 @@ std::vector<double> const& RichardsMechanicsLocalAssembler<
pressure_size> const>(local_x.data() + pressure_index,
pressure_size);
+ auto const material_id =
+ _process_data.flow_material->getMaterialID(_element.getID());
+
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -640,10 +672,11 @@ std::vector<double> const& RichardsMechanicsLocalAssembler<
NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
auto const temperature = _process_data.temperature(t, x_position)[0];
- auto const mu = _process_data.flow_material->getFluidViscosity(
- -p_cap_ip, temperature);
- double const K_over_mu =
- _process_data.intrinsic_permeability(t, x_position)[0] / mu;
+ GlobalDimMatrixType const K_over_mu =
+ intrinsicPermeability<DisplacementDim>(
+ t, x_position, material_id, *_process_data.flow_material) /
+ _process_data.flow_material->getFluidViscosity(-p_cap_ip,
+ temperature);
auto const rho_LR = _process_data.flow_material->getFluidDensity(
-p_cap_ip, temperature);
auto const& b = _process_data.specific_body_force;
@@ -651,7 +684,7 @@ std::vector<double> const& RichardsMechanicsLocalAssembler<
// Compute the velocity
auto const& dNdx_p = _ip_data[ip].dNdx_p;
cache_matrix.col(ip).noalias() =
- -K_over_mu * dNdx_p * p_L - K_over_mu * rho_LR * b;
+ -K_over_mu * dNdx_p * p_L - rho_LR * K_over_mu * b;
}
return cache;
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
index 756705788fb24e340edd595a79956d2f6afb7838..1632ab1a4831e243818819a2166d166054c87cd1 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -49,6 +49,9 @@ public:
using ShapeMatricesTypePressure =
ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>;
+ using GlobalDimMatrixType =
+ typename ShapeMatricesTypePressure::GlobalDimMatrixType;
+
static int const KelvinVectorSize =
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value;
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
index 986309cdf133f608abbb02624a89a8cf2da8687c..6b22b03b1de2a6dc1845f7407a3df33e7a148327 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
@@ -42,7 +42,6 @@ struct RichardsMechanicsProcessData
std::unique_ptr<
MaterialLib::Solids::MechanicsBase<DisplacementDim>>>&&
solid_materials_,
- Parameter<double> const& intrinsic_permeability_,
Parameter<double> const& fluid_bulk_modulus_,
Parameter<double> const& biot_coefficient_,
Parameter<double> const& solid_density_,
@@ -53,7 +52,6 @@ struct RichardsMechanicsProcessData
: material_ids(material_ids_),
flow_material{std::move(flow_material_)},
solid_materials{std::move(solid_materials_)},
- intrinsic_permeability(intrinsic_permeability_),
fluid_bulk_modulus(fluid_bulk_modulus_),
biot_coefficient(biot_coefficient_),
solid_density(solid_density_),
@@ -83,8 +81,6 @@ struct RichardsMechanicsProcessData
int,
std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
solid_materials;
- /// Permeability of the solid. A scalar quantity, Parameter<double>.
- Parameter<double> const& intrinsic_permeability;
/// Fluid's bulk modulus. A scalar quantity, Parameter<double>.
Parameter<double> const& fluid_bulk_modulus;
/// Biot coefficient. A scalar quantity, Parameter<double>.
diff --git a/Tests/Data/RichardsMechanics/RichardsFlow_2d_quasinewton.prj b/Tests/Data/RichardsMechanics/RichardsFlow_2d_quasinewton.prj
index 39d5eb4bf1dcfed4dbd4dbb220895ce9931b7526..e92ab00575c2406c49e3bda627ec3f20dba12e06 100644
--- a/Tests/Data/RichardsMechanics/RichardsFlow_2d_quasinewton.prj
+++ b/Tests/Data/RichardsMechanics/RichardsFlow_2d_quasinewton.prj
@@ -18,7 +18,6 @@
<youngs_modulus>E</youngs_modulus>
<poissons_ratio>nu</poissons_ratio>
</constitutive_relation>
- <intrinsic_permeability>k</intrinsic_permeability>
<solid_bulk_modulus>K_SR</solid_bulk_modulus>
<fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
<biot_coefficient>alpha</biot_coefficient>
diff --git a/Tests/Data/RichardsMechanics/RichardsFlow_2d_small.prj b/Tests/Data/RichardsMechanics/RichardsFlow_2d_small.prj
index 76f932e81a12203bf0d47fbab9502ff7bce69b97..3faa5c9cdad4d8b7015bffebeb6e2e68a3ea1d76 100644
--- a/Tests/Data/RichardsMechanics/RichardsFlow_2d_small.prj
+++ b/Tests/Data/RichardsMechanics/RichardsFlow_2d_small.prj
@@ -13,7 +13,6 @@
<youngs_modulus>E</youngs_modulus>
<poissons_ratio>nu</poissons_ratio>
</constitutive_relation>
- <intrinsic_permeability>k</intrinsic_permeability>
<solid_bulk_modulus>K_SR</solid_bulk_modulus>
<fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
<biot_coefficient>alpha</biot_coefficient>
diff --git a/Tests/Data/RichardsMechanics/confined_compression_fully_saturated.prj b/Tests/Data/RichardsMechanics/confined_compression_fully_saturated.prj
index 803d2011607593dc89e1cd7b91bcc8f5da92cb93..829d6feb2a7b83a9f8540fd57ba8a16c22ca2208 100644
--- a/Tests/Data/RichardsMechanics/confined_compression_fully_saturated.prj
+++ b/Tests/Data/RichardsMechanics/confined_compression_fully_saturated.prj
@@ -20,7 +20,6 @@
<youngs_modulus>E</youngs_modulus>
<poissons_ratio>nu</poissons_ratio>
</constitutive_relation>
- <intrinsic_permeability>k</intrinsic_permeability>
<solid_bulk_modulus>K_SR</solid_bulk_modulus>
<fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
<biot_coefficient>alpha</biot_coefficient>
diff --git a/Tests/Data/RichardsMechanics/flow_fully_saturated.prj b/Tests/Data/RichardsMechanics/flow_fully_saturated.prj
index 0de1a5157f567c7b9fa23cca84d7680cebae2d15..f919f42a3c8285910d92a6ce00bbaa94da99b651 100644
--- a/Tests/Data/RichardsMechanics/flow_fully_saturated.prj
+++ b/Tests/Data/RichardsMechanics/flow_fully_saturated.prj
@@ -13,7 +13,6 @@
<youngs_modulus>E</youngs_modulus>
<poissons_ratio>nu</poissons_ratio>
</constitutive_relation>
- <intrinsic_permeability>k</intrinsic_permeability>
<solid_bulk_modulus>K_SR</solid_bulk_modulus>
<fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
<biot_coefficient>alpha</biot_coefficient>
diff --git a/Tests/Data/RichardsMechanics/gravity.prj b/Tests/Data/RichardsMechanics/gravity.prj
index 1871179477170087e99626a8e0938386b3b7ece5..3e492acb986dc113fd3eac96603874453a9d56b2 100644
--- a/Tests/Data/RichardsMechanics/gravity.prj
+++ b/Tests/Data/RichardsMechanics/gravity.prj
@@ -20,7 +20,6 @@
<youngs_modulus>E</youngs_modulus>
<poissons_ratio>nu</poissons_ratio>
</constitutive_relation>
- <intrinsic_permeability>k</intrinsic_permeability>
<solid_bulk_modulus>K_SR</solid_bulk_modulus>
<fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
<biot_coefficient>alpha</biot_coefficient>
@@ -43,7 +42,7 @@
<porous_medium>
<porous_medium id="0">
<permeability>
- <permeability_tensor_entries>kappa1</permeability_tensor_entries>
+ <permeability_tensor_entries>k</permeability_tensor_entries>
<type>Constant</type>
</permeability>
<porosity>
@@ -200,11 +199,6 @@
<type>Constant</type>
<value>293.15</value>
</parameter>
- <parameter>
- <name>kappa1</name>
- <type>Constant</type>
- <values>4.46e-13</values>
- </parameter>
<parameter>
<name>displacement0</name>
<type>Constant</type>
diff --git a/Tests/Data/RichardsMechanics/mechanics_linear.prj b/Tests/Data/RichardsMechanics/mechanics_linear.prj
index e38e3d1fe87787aaaa76bd086690432b07866301..74dfa7758729a16dedec7a434c589f26f3bcb6f3 100644
--- a/Tests/Data/RichardsMechanics/mechanics_linear.prj
+++ b/Tests/Data/RichardsMechanics/mechanics_linear.prj
@@ -26,7 +26,6 @@
<youngs_modulus>E</youngs_modulus>
<poissons_ratio>nu</poissons_ratio>
</constitutive_relation>
- <intrinsic_permeability>k</intrinsic_permeability>
<solid_bulk_modulus>K_SR</solid_bulk_modulus>
<fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
<biot_coefficient>alpha</biot_coefficient>