Commit 7acf4160 authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov
Browse files

[MPL] Replace volumetric strain rate w/ increments

Add extra comment on the misuse of the volumetric strain
for the mechanical volumetric strain.
parent 79b1dcd0
......@@ -57,8 +57,11 @@ PropertyDataType PorosityFromMassBalance::value(
->property(PropertyType::biot_coefficient)
.template value<double>(variable_array, pos, t, dt);
double const e_dot = std::get<double>(
variable_array[static_cast<int>(Variable::volumetric_strain_rate)]);
double const e = std::get<double>(
variable_array[static_cast<int>(Variable::volumetric_strain)]);
double const e_prev = std::get<double>(
variable_array_prev[static_cast<int>(Variable::volumetric_strain)]);
double const delta_e = e - e_prev;
double const p_eff = std::get<double>(
variable_array[static_cast<int>(Variable::effective_pore_pressure)]);
......@@ -70,7 +73,7 @@ PropertyDataType PorosityFromMassBalance::value(
double const phi_prev = std::get<double>(
variable_array_prev[static_cast<int>(Variable::porosity)]);
double const w = dt * e_dot + delta_p_eff * beta_SR;
double const w = delta_e + delta_p_eff * beta_SR;
return std::clamp((phi_prev + alpha_b * w) / (1 + w), phi_min_, phi_max_);
}
......
......@@ -47,8 +47,11 @@ PropertyDataType TransportPorosityFromMassBalance::value(
->property(PropertyType::biot_coefficient)
.template value<double>(variable_array, pos, t, dt);
double const e_dot = std::get<double>(
variable_array[static_cast<int>(Variable::volumetric_strain_rate)]);
double const e = std::get<double>(
variable_array[static_cast<int>(Variable::volumetric_strain)]);
double const e_prev = std::get<double>(
variable_array_prev[static_cast<int>(Variable::volumetric_strain)]);
double const delta_e = e - e_prev;
double const p_eff = std::get<double>(
variable_array[static_cast<int>(Variable::effective_pore_pressure)]);
......@@ -63,7 +66,7 @@ PropertyDataType TransportPorosityFromMassBalance::value(
double const phi_tr_prev = std::get<double>(
variable_array_prev[static_cast<int>(Variable::transport_porosity)]);
double const w = dt * e_dot + delta_p_eff * beta_SR;
double const w = delta_e + delta_p_eff * beta_SR;
return std::clamp(phi_tr_prev + (alpha_b - phi) * w, phi_min_, phi_max_);
}
......
......@@ -56,7 +56,7 @@ enum class Variable : int
stress,
temperature,
transport_porosity,
volumetric_strain_rate,
volumetric_strain,
number_of_variables
};
......
......@@ -387,9 +387,10 @@ void RichardsMechanicsLocalAssembler<
-chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
// Set volumetric strain rate for the general case without swelling.
double const div_u_dot = identity2.transpose() * B * u_dot;
variables[static_cast<int>(MPL::Variable::volumetric_strain_rate)]
.emplace<double>(div_u_dot);
variables[static_cast<int>(MPL::Variable::volumetric_strain)]
.emplace<double>(identity2.transpose() * B * u);
variables_prev[static_cast<int>(MPL::Variable::volumetric_strain)]
.emplace<double>(identity2.transpose() * B * (u - u_dot * dt));
auto& phi = _ip_data[ip].porosity;
{ // Porosity update
......@@ -424,11 +425,14 @@ void RichardsMechanicsLocalAssembler<
auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
t, x_position, dt, temperature);
variables[static_cast<int>(
MPL::Variable::volumetric_strain_rate)]
.emplace<double>(div_u_dot + identity2.transpose() *
C_el.inverse() *
sigma_sw_dot);
// !!! 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 (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
......@@ -715,9 +719,10 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
-chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
// Set volumetric strain rate for the general case without swelling.
double const div_u_dot = identity2.transpose() * B * u_dot;
variables[static_cast<int>(MPL::Variable::volumetric_strain_rate)]
.emplace<double>(div_u_dot);
variables[static_cast<int>(MPL::Variable::volumetric_strain)]
.emplace<double>(identity2.transpose() * B * u);
variables_prev[static_cast<int>(MPL::Variable::volumetric_strain)]
.emplace<double>(identity2.transpose() * B * (u - u_dot * dt));
auto& phi = _ip_data[ip].porosity;
{ // Porosity update
......@@ -750,11 +755,14 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables, variables_prev, x_position, t, dt));
sigma_sw += sigma_sw_dot * dt;
variables[static_cast<int>(
MPL::Variable::volumetric_strain_rate)]
.emplace<double>(div_u_dot + identity2.transpose() *
C_el.inverse() *
sigma_sw_dot);
// !!! 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 (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
......@@ -1441,9 +1449,10 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
-chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
// Set volumetric strain rate for the general case without swelling.
double const div_u_dot = identity2.transpose() * B * u_dot;
variables[static_cast<int>(MPL::Variable::volumetric_strain_rate)]
.emplace<double>(div_u_dot);
variables[static_cast<int>(MPL::Variable::volumetric_strain)]
.emplace<double>(identity2.transpose() * B * u);
variables_prev[static_cast<int>(MPL::Variable::volumetric_strain)]
.emplace<double>(identity2.transpose() * B * (u - u_dot * dt));
auto& phi = _ip_data[ip].porosity;
{ // Porosity update
......@@ -1475,11 +1484,14 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables, variables_prev, x_position, t, dt));
sigma_sw += sigma_sw_dot * dt;
variables[static_cast<int>(
MPL::Variable::volumetric_strain_rate)]
.emplace<double>(div_u_dot + identity2.transpose() *
C_el.inverse() *
sigma_sw_dot);
// !!! 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 (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
......
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