From c9de7b2bd0cee7aaf6e62edb646c99c751ddcd50 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Wed, 26 Feb 2020 17:16:41 +0100
Subject: [PATCH] [PL/RM] Update porosity computation.
Now includes the Bishop's effective stress.
---
.../RichardsMechanicsFEM-impl.h | 65 ++++++++++++-------
1 file changed, 41 insertions(+), 24 deletions(-)
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 217a13d2247..c1368c241a0 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -251,11 +251,12 @@ void RichardsMechanicsLocalAssembler<
double p_cap_ip;
NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
+ double p_cap_dot_ip;
+ NumLib::shapeFunctionInterpolate(-p_L_dot, N_p, p_cap_dot_ip);
+
variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
p_cap_ip;
variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
- variables[static_cast<int>(MPL::Variable::phase_pressure_rate)] =
- N_p.dot(p_L_dot);
auto const temperature =
medium->property(MPL::PropertyType::reference_temperature)
@@ -275,15 +276,6 @@ void RichardsMechanicsLocalAssembler<
liquid_phase.property(MPL::PropertyType::bulk_modulus)
.template value<double>(variables, x_position, t, dt);
- // Use previous time step porosity for porosity update, ...
- variables[static_cast<int>(MPL::Variable::porosity)] =
- _ip_data[ip].porosity_prev;
- auto& porosity = _ip_data[ip].porosity;
- porosity = solid_phase.property(MPL::PropertyType::porosity)
- .template value<double>(variables, x_position, t, dt);
- // ... then use new porosity.
- variables[static_cast<int>(MPL::Variable::porosity)] = porosity;
-
auto const rho_LR =
liquid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
@@ -302,9 +294,30 @@ void RichardsMechanicsLocalAssembler<
MPL::Variable::capillary_pressure,
x_position, t, dt);
- auto const chi_S_L =
- medium->property(MPL::PropertyType::bishops_effective_stress)
+ auto const chi = [&](double const S_L) {
+ MPL::VariableArray variables;
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
+ return medium->property(MPL::PropertyType::bishops_effective_stress)
.template value<double>(variables, x_position, t, dt);
+ };
+ double const chi_S_L = chi(S_L);
+ double const chi_S_L_prev = chi(S_L_prev);
+
+ variables[static_cast<int>(
+ MPL::Variable::effective_pore_pressure_rate)] =
+ (chi_S_L * (-p_cap_ip) -
+ chi_S_L_prev * (-p_cap_ip + p_cap_dot_ip * dt)) /
+ dt;
+
+ // Use previous time step porosity for porosity update, ...
+ variables[static_cast<int>(MPL::Variable::porosity)] =
+ _ip_data[ip].porosity_prev;
+ auto& porosity = _ip_data[ip].porosity;
+ porosity = solid_phase.property(MPL::PropertyType::porosity)
+ .template value<double>(variables, x_position, t, dt);
+ // ... then use new porosity.
+ variables[static_cast<int>(MPL::Variable::porosity)] = porosity;
+
double const k_rel =
medium->property(MPL::PropertyType::relative_permeability)
@@ -516,8 +529,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
p_cap_ip;
variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
- variables[static_cast<int>(MPL::Variable::phase_pressure_rate)] =
- N_p.dot(p_L_dot);
variables[static_cast<int>(MPL::Variable::volumetric_strain_rate)]
.emplace<double>(identity2.transpose() * B * u_dot);
auto const temperature =
@@ -544,15 +555,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
liquid_phase.property(MPL::PropertyType::bulk_modulus)
.template value<double>(variables, x_position, t, dt);
- // Use previous time step porosity for porosity update, ...
- variables[static_cast<int>(MPL::Variable::porosity)] =
- _ip_data[ip].porosity_prev;
- auto& porosity = _ip_data[ip].porosity;
- porosity = solid_phase.property(MPL::PropertyType::porosity)
- .template value<double>(variables, x_position, t, dt);
- // ... then use new porosity.
- variables[static_cast<int>(MPL::Variable::porosity)] = porosity;
-
auto const rho_LR =
liquid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
@@ -596,6 +598,21 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu;
+ variables[static_cast<int>(
+ MPL::Variable::effective_pore_pressure_rate)] =
+ (chi_S_L * (-p_cap_ip) -
+ chi_S_L_prev * (-p_cap_ip + p_cap_dot_ip * dt)) /
+ dt;
+
+ // Use previous time step porosity for porosity update, ...
+ variables[static_cast<int>(MPL::Variable::porosity)] =
+ _ip_data[ip].porosity_prev;
+ auto& porosity = _ip_data[ip].porosity;
+ porosity = solid_phase.property(MPL::PropertyType::porosity)
+ .template value<double>(variables, x_position, t, dt);
+ // ... then use new porosity.
+ variables[static_cast<int>(MPL::Variable::porosity)] = porosity;
+
sigma_sw = sigma_sw_prev;
if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
{
--
GitLab