Commit 4e230d0c authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov
Browse files

Merge branch 'RMMicroporosityFixes' into 'master'

[PL/RM] Reuse swelling stress update in secondary variables update.

See merge request ogs/ogs!3585
parents 00a14b84 f40334c5
......@@ -66,7 +66,7 @@ MicroPorosityStateSpace<DisplacementDim> computeMicroPorosity(
double const rho_LR_m, // for simplification equal to rho_LR_M
double const mu_LR,
MicroPorosityParameters const& micro_porosity_parameters,
double const alpha_B, double const phi_M, double const p_L,
double const alpha_B, double const phi, double const p_L,
double const p_L_m_prev,
MaterialPropertyLib::VariableArray const& /*variables_prev*/,
double const S_L_m_prev, double const phi_m_prev,
......@@ -135,14 +135,14 @@ MicroPorosityStateSpace<DisplacementDim> computeMicroPorosity(
swelling_stress_rate.template value<Eigen::Matrix3d>(
variables, variables_prev, pos, t, dt));
residual[i_phi_m] = delta_phi_m - (alpha_B - phi_M) * delta_e_sw;
residual[i_phi_m] = delta_phi_m - (alpha_B - phi) * delta_e_sw;
residual[i_e_sw] = delta_e_sw + I_2_C_el_inverse.dot(delta_sigma_sw);
residual.template segment<kelvin_vector_size>(i_sigma_sw).noalias() =
delta_sigma_sw - sigma_sw_dot * dt;
residual[i_p_L_m] =
rho_LR_m *
(phi_m * delta_S_L_m - (alpha_B - phi_M) * S_L_m * delta_e_sw) +
(phi_m * delta_S_L_m - (alpha_B - phi) * S_L_m * delta_e_sw) +
phi_m * S_L_m * rho_LR_m * delta_e_sw -
micro_porosity_parameters.mass_exchange_coefficient / mu_LR *
(p_L - p_L_m) * dt;
......@@ -178,7 +178,7 @@ MicroPorosityStateSpace<DisplacementDim> computeMicroPorosity(
variables, variables_prev, MPL::Variable::liquid_saturation,
pos, t, dt));
jacobian(i_phi_m, i_e_sw) = -(alpha_B - phi_M);
jacobian(i_phi_m, i_e_sw) = -(alpha_B - phi);
jacobian.template block<1, kelvin_vector_size>(i_e_sw, i_sigma_sw) =
I_2_C_el_inverse.transpose();
......@@ -189,11 +189,11 @@ MicroPorosityStateSpace<DisplacementDim> computeMicroPorosity(
jacobian(i_p_L_m, i_phi_m) =
rho_LR_m * (delta_S_L_m + S_L_m * delta_e_sw);
jacobian(i_p_L_m, i_e_sw) = -rho_LR_m * S_L_m * (alpha_B - phi_M - phi_m);
jacobian(i_p_L_m, i_e_sw) = -rho_LR_m * S_L_m * (alpha_B - phi - phi_m);
jacobian(i_p_L_m, i_p_L_m) =
alpha_bar / mu_LR * dt -
rho_LR_m * (phi_m - (alpha_B - phi_M - phi_m) * delta_e_sw) *
rho_LR_m * (phi_m - (alpha_B - phi - phi_m) * delta_e_sw) *
dS_L_m_dp_cap_m;
};
......
......@@ -27,6 +27,94 @@ namespace ProcessLib
{
namespace RichardsMechanics
{
template <int DisplacementDim, typename IPData>
void updateSwellingStressAndVolumetricStrain(
IPData& ip_data, MaterialPropertyLib::Medium const& medium,
MaterialPropertyLib::Phase const& solid_phase,
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> const& C_el,
double const rho_LR, double const mu,
std::optional<MicroPorosityParameters> micro_porosity_parameters,
double const alpha, double const phi, double const p_cap_ip,
MPL::VariableArray& variables, MPL::VariableArray& variables_prev,
ParameterLib::SpatialPosition const& x_position, double const t,
double const dt)
{
auto const& identity2 = MathLib::KelvinVector::Invariants<
MathLib::KelvinVector::kelvin_vector_dimensions(
DisplacementDim)>::identity2;
auto& sigma_sw = ip_data.sigma_sw;
auto const& sigma_sw_prev = ip_data.sigma_sw_prev;
if (!medium.hasProperty(MPL::PropertyType::saturation_micro))
{
// If there is swelling, compute it. Update volumetric strain rate,
// s.t. it corresponds to the mechanical part only.
sigma_sw = sigma_sw_prev;
if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
{
using DimMatrix = Eigen::Matrix<double, 3, 3>;
auto const sigma_sw_dot =
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
solid_phase
.property(MPL::PropertyType::swelling_stress_rate)
.template value<DimMatrix>(variables, variables_prev,
x_position, t, dt));
sigma_sw += sigma_sw_dot * dt;
// !!! Misusing volumetric strain for mechanical volumetric
// strain just to update the transport porosity !!!
std::get<double>(variables[static_cast<int>(
MPL::Variable::volumetric_strain)]) +=
identity2.transpose() * C_el.inverse() * sigma_sw;
std::get<double>(variables_prev[static_cast<int>(
MPL::Variable::volumetric_strain)]) +=
identity2.transpose() * C_el.inverse() * sigma_sw_prev;
}
}
// 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.transport_porosity_prev;
double const phi_prev = ip_data.porosity_prev;
double const phi_m_prev = phi_prev - phi_M_prev;
double const p_L_m_prev = ip_data.liquid_pressure_m_prev;
auto const S_L_m_prev = ip_data.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,
*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.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.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.saturation_m =
medium.property(MPL::PropertyType::saturation_micro)
.template value<double>(variables, x_position, t, dt);
}
sigma_sw = sigma_sw_prev + delta_sigma_sw;
}
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
......@@ -837,86 +925,15 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
alpha, phi, _element.getID(), ip);
}
// 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))
{
// If there is swelling, compute it. Update volumetric strain rate,
// s.t. it corresponds to the mechanical part only.
sigma_sw = sigma_sw_prev;
if (solid_phase.hasProperty(
MPL::PropertyType::swelling_stress_rate))
{
using DimMatrix = Eigen::Matrix<double, 3, 3>;
auto const sigma_sw_dot =
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
solid_phase
.property(MPL::PropertyType::swelling_stress_rate)
.template value<DimMatrix>(
variables, variables_prev, x_position, t, dt));
sigma_sw += sigma_sw_dot * dt;
// !!! Misusing volumetric strain for mechanical volumetric
// strain just to update the transport porosity !!!
std::get<double>(variables[static_cast<int>(
MPL::Variable::volumetric_strain)]) +=
identity2.transpose() * C_el.inverse() * sigma_sw;
std::get<double>(variables_prev[static_cast<int>(
MPL::Variable::volumetric_strain)]) +=
identity2.transpose() * C_el.inverse() * sigma_sw_prev;
}
}
auto const mu =
liquid_phase.property(MPL::PropertyType::viscosity)
.template value<double>(variables, x_position, t, dt);
// 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;
}
// Swelling and possibly volumetric strain rate update.
updateSwellingStressAndVolumetricStrain<DisplacementDim>(
_ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu,
_process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
variables, variables_prev, x_position, t, dt);
if (medium->hasProperty(MPL::PropertyType::transport_porosity))
{
......@@ -966,6 +983,8 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
//
// displacement equation, displacement part
//
auto& sigma_sw = _ip_data[ip].sigma_sw;
eps_m.noalias() =
solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
? eps + C_el.inverse() * sigma_sw
......@@ -1687,55 +1706,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables[static_cast<int>(MPL::Variable::porosity)] = phi;
}
// 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 there is swelling, compute it. Update volumetric strain rate,
// s.t. it corresponds to the mechanical part only.
sigma_sw = sigma_sw_prev;
if (solid_phase.hasProperty(
MPL::PropertyType::swelling_stress_rate))
{
using DimMatrix = Eigen::Matrix<double, 3, 3>;
auto const sigma_sw_dot =
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
solid_phase
.property(MPL::PropertyType::swelling_stress_rate)
.template value<DimMatrix>(
variables, variables_prev, x_position, t, dt));
sigma_sw += sigma_sw_dot * dt;
// !!! Misusing volumetric strain for mechanical volumetric
// strain just to update the transport porosity !!!
std::get<double>(variables[static_cast<int>(
MPL::Variable::volumetric_strain)]) +=
identity2.transpose() * C_el.inverse() * sigma_sw;
std::get<double>(variables_prev[static_cast<int>(
MPL::Variable::volumetric_strain)]) +=
identity2.transpose() * C_el.inverse() * sigma_sw_prev;
}
if (medium->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 =
medium->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;
}
}
auto const mu =
liquid_phase.property(MPL::PropertyType::viscosity)
.template value<double>(variables, x_position, t, dt);
......@@ -1743,6 +1713,31 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
liquid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
// Swelling and possibly volumetric strain rate update.
updateSwellingStressAndVolumetricStrain<DisplacementDim>(
_ip_data[ip], *medium, solid_phase, C_el, rho_LR, mu,
_process_data.micro_porosity_parameters, alpha, phi, p_cap_ip,
variables, variables_prev, x_position, t, dt);
if (medium->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 =
medium->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;
}
// Set mechanical variables for the intrinsic permeability model
// For stress dependent permeability.
{
......@@ -1780,6 +1775,8 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
.template value<double>(variables, x_position, t, dt);
_ip_data[ip].dry_density_solid = (1 - phi) * rho_SR;
auto& sigma_sw = _ip_data[ip].sigma_sw;
eps_m.noalias() =
solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
? eps + C_el.inverse() * sigma_sw
......
......@@ -334,7 +334,7 @@
<vtkdiff>
<regex>double_porosity_swelling_.*.vtu</regex>
<field>sigma</field>
<absolute_tolerance>1e-8</absolute_tolerance>
<absolute_tolerance>2e-8</absolute_tolerance>
<relative_tolerance>0</relative_tolerance>
</vtkdiff>
<vtkdiff>
......@@ -365,7 +365,7 @@
<regex>double_porosity_swelling_.*.vtu</regex>
<field>dry_density_solid</field>
<absolute_tolerance>2e-12</absolute_tolerance>
<relative_tolerance>1e-15</relative_tolerance>
<relative_tolerance>2e-15</relative_tolerance>
</vtkdiff>
</test_definition>
</OpenGeoSysProject>
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment