diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 7a1dce59408f9dfd53d83e03349c72a443062ca7..98f57b37cdc1f97058b1956617413004a44ca0d9 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -13,6 +13,7 @@
#include <cassert>
+#include "ComputeMicroPorosity.h"
#include "IntegrationPointData.h"
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
@@ -231,7 +232,16 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
medium->property(MPL::PropertyType::saturation)
.template value<double>(variables, x_position, t, dt);
- _ip_data[ip].saturation_m_prev = _ip_data[ip].saturation_prev;
+ if (medium->hasProperty(MPL::PropertyType::saturation_micro))
+ {
+ MPL::VariableArray vars;
+ vars[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ p_cap_ip;
+ double const S_L_m =
+ medium->property(MPL::PropertyType::saturation_micro)
+ .template value<double>(vars, x_position, t, dt);
+ _ip_data[ip].saturation_m_prev = S_L_m;
+ }
// Set eps_m_prev from potentially non-zero eps and sigma_sw from
// restart.
@@ -805,8 +815,9 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
// Swelling and possibly volumetric strain rate update.
auto& sigma_sw = _ip_data[ip].sigma_sw;
+ auto const& sigma_sw_prev = _ip_data[ip].sigma_sw_prev;
+ if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
{
- auto const& sigma_sw_prev = _ip_data[ip].sigma_sw_prev;
// If there is swelling, compute it. Update volumetric strain rate,
// s.t. it corresponds to the mechanical part only.
@@ -832,33 +843,79 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
MPL::Variable::volumetric_strain)]) +=
identity2.transpose() * C_el.inverse() * sigma_sw_prev;
}
+ }
- if (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
- {
- variables_prev[static_cast<int>(
- MPL::Variable::transport_porosity)] =
- _ip_data[ip].transport_porosity_prev;
+ auto const mu =
+ liquid_phase.property(MPL::PropertyType::viscosity)
+ .template value<double>(variables, x_position, t, dt);
- _ip_data[ip].transport_porosity =
- solid_phase.property(MPL::PropertyType::transport_porosity)
- .template value<double>(variables, variables_prev,
- x_position, t, dt);
- variables[static_cast<int>(MPL::Variable::transport_porosity)] =
- _ip_data[ip].transport_porosity;
- }
- else
- {
- variables[static_cast<int>(MPL::Variable::transport_porosity)] =
- phi;
+ // TODO (naumov) saturation_micro must be always defined together with
+ // the micro_porosity_parameters.
+ if (medium->hasProperty(MPL::PropertyType::saturation_micro))
+ {
+ double const phi_M_prev = _ip_data[ip].transport_porosity_prev;
+ double const phi_prev = _ip_data[ip].porosity_prev;
+ double const phi_m_prev = phi_prev - phi_M_prev;
+ double const p_L_m_prev = _ip_data[ip].liquid_pressure_m_prev;
+
+ auto const S_L_m_prev = _ip_data[ip].saturation_m_prev;
+
+ auto const [delta_phi_m, delta_e_sw, delta_p_L_m, delta_sigma_sw] =
+ computeMicroPorosity<DisplacementDim>(
+ identity2.transpose() * C_el.inverse(), rho_LR, mu,
+ *_process_data.micro_porosity_parameters, alpha, phi,
+ -p_cap_ip, p_L_m_prev, variables_prev, S_L_m_prev,
+ phi_m_prev, x_position, t, dt,
+ medium->property(MPL::PropertyType::saturation_micro),
+ solid_phase.property(
+ MPL::PropertyType::swelling_stress_rate));
+
+ auto& phi_M = _ip_data[ip].transport_porosity;
+ phi_M = phi - (phi_m_prev + delta_phi_m);
+ variables_prev[static_cast<int>(
+ MPL::Variable::transport_porosity)] = phi_M_prev;
+ variables[static_cast<int>(MPL::Variable::transport_porosity)] =
+ phi_M;
+
+ auto& p_L_m = _ip_data[ip].liquid_pressure_m;
+ p_L_m = p_L_m_prev + delta_p_L_m;
+ { // Update micro saturation.
+ MPL::VariableArray variables_prev;
+ variables_prev[static_cast<int>(
+ MPL::Variable::capillary_pressure)] = -p_L_m_prev;
+ MPL::VariableArray variables;
+ variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ -p_L_m;
+
+ _ip_data[ip].saturation_m =
+ medium->property(MPL::PropertyType::saturation_micro)
+ .template value<double>(variables, x_position, t, dt);
}
+ sigma_sw = sigma_sw_prev + delta_sigma_sw;
+ }
+
+ if (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
+ {
+ variables_prev[static_cast<int>(
+ MPL::Variable::transport_porosity)] =
+ _ip_data[ip].transport_porosity_prev;
+
+ _ip_data[ip].transport_porosity =
+ solid_phase.property(MPL::PropertyType::transport_porosity)
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::transport_porosity)] =
+ _ip_data[ip].transport_porosity;
+ }
+ else
+ {
+ variables[static_cast<int>(MPL::Variable::transport_porosity)] =
+ phi;
}
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);
// Set mechanical variables for the intrinsic permeability model
// For stress dependent permeability.
@@ -938,7 +995,17 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
.noalias() +=
N_u_op.transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
- if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
+ if (medium->hasProperty(MPL::PropertyType::saturation_micro))
+ {
+ // For the swelling stress with double structure model the
+ // corresponding Jacobian u-p entry would be:
+ // -B.transpose() *
+ // dsigma_sw_dS_L_m* dS_L_m_dp_cap_m*(p_L_m - p_L_m_prev) /
+ // p_cap_dot_ip / dt* N_p* w;
+ // but it does not improve convergence and sometimes worsen it.
+ }
+ else if (solid_phase.hasProperty(
+ MPL::PropertyType::swelling_stress_rate))
{
using DimMatrix = Eigen::Matrix<double, 3, 3>;
auto const dsigma_sw_dS_L =
@@ -1040,6 +1107,31 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
+
+ if (medium->hasProperty(MPL::PropertyType::saturation_micro))
+ {
+ double const alpha_bar = _process_data.micro_porosity_parameters
+ ->mass_exchange_coefficient;
+ auto const p_L_m = _ip_data[ip].liquid_pressure_m;
+ local_rhs.template segment<pressure_size>(pressure_index)
+ .noalias() -=
+ N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w;
+
+ local_Jac
+ .template block<pressure_size, pressure_size>(pressure_index,
+ pressure_index)
+ .noalias() += N_p.transpose() * alpha_bar / mu * N_p * w;
+ if (p_cap_dot_ip != 0)
+ {
+ double const p_L_m_prev = _ip_data[ip].liquid_pressure_m_prev;
+ local_Jac
+ .template block<pressure_size, pressure_size>(
+ pressure_index, pressure_index)
+ .noalias() += N_p.transpose() * alpha_bar / mu *
+ (p_L_m - p_L_m_prev) / (dt * p_cap_dot_ip) *
+ N_p * w;
+ }
+ }
}
if (_process_data.apply_mass_lumping)